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ABSTRACT 


A  solution  methodology  for  the  numerical  simulation  of  turbulent  reacting  flows  in 
geometries  described  by  general  curvilinear  coordinates  has  been  developed.  The  method 
is  based  on  the  solution  of  Favre  averaged  variables.  A  k-e  turbulence  model  with  wall 

functions  has  been  used  to  close  the  governing  equations.  Models  for  turbulent  diffusion 
and  premixed  flames  have  been  incorporated  into  the  solution  method.  The  numerical 
solution  of  the  governing  equations  is  accomplished  by  solving  for  the  Cartesian  velocity 
components  on  a  non-staggered  grid.  An  evaluation  of  the  solution  methodology  is 
presented. 
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NOMENCLATURE 


ail 

ai2 

a2i 

a22 

Ap,Ae,  Aw 
An»  As 
C 

Cebu 

Ci 

C2 

Cp 

D 

E 

Ea 

f 

fx 

fy 

g 

h 

h0 

J 

k 

kf 


Transformation  metric;  an  =  y^ 

Transformation  metric;  ai2  =  -x^  _ 

Transformation  metric;  a2i  =  -y% 

Transformation  metric;  &22  =  x^. 

Coefficients  in  the  discrete  equation  for  a  given  variable. 

Molar  concentration  (mol/cnP). 

Constant  in  eddy  break  up  model  of  turbulent  combustion,  3.0. 
Turbulence  constant,  0.09. 

Turbulence  constant,  1.44. 

Turbulence  constant,  1.92. 

Constant  pressure  specific  heat  (J/kgK). 

Diameter. 

Constant  in  log-law  relation,  9.025. 

Activation  energy  (J/kmol). 

Mixture  fraction. 

Interpolation  factor  for  ^  coordinate. 

Interpolation  factor  for  T)  coordinate. 

Concentration  fluctuation. 

Enthalpy  (J/kg). 

Stagnation  enthalpy  (J/kg). 

Jacobian  of  the  coordinate  transformation. 

Turbulent  kinetic  energy  (m2/s2). 

Reaction  rate  constant 


vrn 


Lt 

nx 

ny 

P 

Pk 

Pr 

qn,qi2 

q2l»q22 

Q 


r,z 


R 

R, 

S 


Sc 

t 

T 

u 

U 

v 

V 

w 

W 


x,y 

Y 


Length  scale  of  turbulence. 

Number  of  cells  in  the  £  direction. 

Number  of  cells  in  the  rj  direction. 

Pressure  (N/m2). 

Production  of  turbulent  kinetic  energy. 

Prandtl  number. 

Metrics  of  the  coordinate  transformation. 

Heat  of  reaction  (J/kg). 

Cartesian  coordinates  used  to  describe  axisymmetric  geometries. 

Gas  constant  for  a  particular  gas.  Also,  the  radius  of  a  pipe. 
Universal  gas  constant,  8314.3  J/kmol  K. 

Designation  given  to  a  source  term  in  the  governing  equations.  Also 
constant  in  Sutherland's  law  expression  for  molecular  viscosity. 
Schmidt  number. 

Time  (s). 

Temperature  (K). 

Cartesian  velocity  component  in  the  x  (z)  direction. 

Contravariant  velocity  component  petpendicular  to  lines  of  constant  £. 
Cartesian  velocity  component  in  the  y  (r)  direction. 

Contravariant  velocity  component  perpendicular  to  lines  of  constant  rj. 
Swirl  velocity  component 
Molecular  weight  (kg/kmol). 

Cartesian  coordinates. 

Mass  fraction. 


Greek  Symbols 


a  Variable  used  in  turbulent  diffusion  flame  model.  Fraction  of  time  spent  in 

a  particular  state. 

8ij  Kronecker  delta  function, 

e  Turbulent  dissipation  rate  (m2/s3). 

<t>  General  scalar  variable.  Also,  the  stoichiometric  air  to  fuel  ratio. 

T  Diffusive  exchange  coefficient 

T|  General  curvilinear  coordinate. 

K  von  Kantian  constant,  0.4. 

(i  Dynamic  molecular  viscosity  (kg/m  s). 

v  Kinematic  viscosity  (m2/s). 

p  Density  (kg/m3), 

x  Stress  acting  on  fluid  (N/m2). 

co  Relaxation  factor, 

cb  Reaction  rate  (kg/m3  s). 

£  General  curvilinear  coordinate. 

VF  Schvab-Zeldovich  variable. 

Subscripts 

ARR  Arrhenius  rate  expression. 

EBU  Eddy  break  up  expression. 

e  Quantity  corresponding  to  the  east  face  of  a  control  volume  ( identical  to  use 

ofi+1/2). 

eff  Effective  value  of  a  quantity  for  a  particular  variable. 


x 


FU  Fuel  variable. 

i  Denotes  a  particular  specie.  Also  the  index  in  £  direction. 

IN  Refers  to  an  inlet  quantity, 

j  Index  in  the  fl  direction. 

L  Denotes  a  laminar  quantity. 

n  Quantity  corresponding  to  the  north  face  of  a  control  volume  ( identical  to 

use  of  j  - 1/2 ). 

nb  Refers  to  neighbor  quantities. 

OX  Oxidant  variable. 

PR  Product  variable. 

s  Quantity  corresponding  to  the  south  face  of  a  control  volume  ( identical  to 

use  of  j  +  1/2 ). 

T  Denotes  turbulence  quantity. 

TR  Denotes  quantity  arising  from  coordinate  transformation, 

w  Quantity  corresponding  to  the  west  face  of  a  control  volume  ( identical  to 

use  of  i  - 1/2 ). 


Superscripts 


Denotes  Reynolds  average  fluctuation.  Also  a  correction  quantity. 
Denotes  Favre  average  fluctuation. 

+  Denotes  a  parameter  non-dimensionalized  with  respect  to  friction  velocity. 

Also  denotes  a  quantity  at  the  upper  limit  of  turbulent  fluctuation. 

Denotes  a  Reynolds  averaged  quantity. 

Also  denotes  a  quantity  at  the  lower  limit  of  turbulent  fluctuation. 
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Denotes  a  Favre  averaged  quantity. 

Denotes  an  estimated  value  of  the  prescribed  quantity. 


CHAPTER  1:  INTRODUCTION 


The  ultimate  goal  in  understanding  combustion  processes  in  practical  devices  is  to  relate 
the  parameters  in  the  control  of  the  designer  to  the  performance  of  the  combustion  system 
(Edelman  and  Harsha,  1978).  Identification  of  the  parameters  which  control  the 
distribution  and  extent  of  heat  release  is  important  in  the  design  of  efficient  combustion 
devices.  It  is  also  desirable  to  know  how  flow  turbulence  effects  the  mixing  and 
combustion  of  fuel  and  air.  The  numerical  simulation  of  turbulent  reacting  flows  using 
computational  fluid  dynamics  is  useful  in  attaining  these  goals. 

A  complete  simulation  of  a  turbulent  reacting  flow  requires  the  solution  of  the  time 
dependent  Navier-Stokes  equations  along  with  the  governing  equations  for  species  and 
energy.  Such  a  simulation  is  referred  to  as  a  direct  numerical  simulation  (DNS)  (Jou  and 
Riley,  1989).  A  direct  numerical  simulation  is  capable  of  resolving  the  temporal  and  spatial 
dependencies  of  the  turbulent  flow  field.  However,  the  ability  to  resolve  the  significant 
time  and  length  scales  is  achieved  at  the  cost  of  computational  effort.  As  an  example,  a 
fully  turbulent  isothermal  channel  flow  calculation  using  a  DNS  required  250  hours  of 
CPU  time  on  a  Cray  X-MP  supercomputer  (Moin  and  Kim,  1982).  At  the  present  time, 
DNS  is  not  a  viable  option  for  the  solution  of  practical  flow  fields  (Pope,  1990). 

In  an  attempt  to  circumvent  the  problem  of  large  computational  times  for  solution,  a 
large  eddy  simulation  (LES)  may  be  performed  (Tafti  and  Vanka,  1990).  The  use  of  LES 
in  simulating  turbulent  reacting  flow  has  been  undertaken  by  some  researchers  (Menon  and 
Jou,  1991).  This  method  appears  to  be  quite  promising  but  computational  times  are  still 
quite  large. 

The  solution  of  time  or  mass  averaged  equations  is  currently  the  most  widely  used 
means  of  simulating  a  turbulent  reacting  flow  (Jou  and  Riley,  1989).  In  solving  a  set  of 
averaged  equations,  a  temporally  and  spatially  resolved  solution  is  replaced  by  a  stochastic 
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description  of  the  flow.  The  dependent  variables  of  the  governing  equations  are  time  or 
mass  averaged  quantities. 

The  solution  of  averaged  equations  is  obtained  only  after  certain  closure  hypotheses 
have  been  invoked.  The  closure  hypotheses  or  turbulence  models  relate  the  mean 
fluctuations  (Reynolds  stresses)  to  the  time  or  mass  averaged  flow  variables.  The  degree 
of  sophistication  of  this  closure  is  somewhat  dependent  on  the  flowfields  which  are  to  be 
studied.  A  descriptive  review  of  the  various  models  which  are  in  use  for  predicting  internal 
flows  is  given  by  Nallasamy,  1987.  For  wall  bounded  flows,  Patel  et.  al.,  1985  give  and 
analysis  of  various  low  Reynolds  number  k-e  turbulence  models.  A  review  of  Reynolds 
stress  models  is  given  by  Speziale,  1991.  A  turbulence  model  of  some  kind  must  be  used 
if  a  solution  of  the  averaged  flow  variables  is  sought. 

In  simulating  turbulent  reacting  flows,  a  combustion  model  must  also  be  used.  The 
objective  of  the  combustion  model  is  to  relate  the  time  mean  reaction  rate  to  the  mean  flow 
variables.  In  general,  different  combustion  models  are  used  depending  on  whether  a 
premixed  or  non-premixed  combustion  process  is  to  be  studied  (Lockwood,  1977). 

For  the  simulation  of  diffusion  or  non-premixed  flames,  it  has  become  common  to 
assume  infinitely  fast  kinetics  and  equal  specie  diffusivities,  whence  the  thermochemical 
state  of  the  mixture  at  any  point  and  time  is  only  dependent  on  the  mixture  fraction  (Bilger, 
1989).  The  averaged  species,  temperature,  and  density  are  determined  either  by  assuming 
a  probability  density  function  (pdf)  (Bilger,  1980)  or  by  solving  an  equation  for  the  pdf 
(Pope,  1990).  As  was  the  case  for  turbulence  models,  different  levels  of  complexity  exist 
among  the  many  models  for  diffusion  flames.  Reviews  by  Spalding,  1976  and  by  Jones 
and  Whitelaw,  1982  present  some  of  the  various  aspects  relevant  to  modelling  a  turbulent 
diffusion  flame. 

Simulation  of  turbulent  premixed  flames  has  proven  to  be  more  difficult  than  the 
simulation  of  diffusion  flames  (Pope,  1990).  A  currently  used  model  for  premixed  flames 
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is  the  eddy  break  up  model  of  Spalding,  1976.  Various  forms  of  this  model  have  been 
proposed  and  are  in  use  (Magnussen  and  Hjertager,  1976;  Lockwood,  1977).  Another 
method  for  modelling  turbulent  premixed  flames  involves  solving  for  a  reaction  progress 
variable.  Details  of  this  method  are  given  by  Bray,  1980. 

Once  the  governing  equations  and  closure  models  have  been  specified,  the  task 
becomes  one  of  solving  a  coupled  set  of  partial  differential  equations.  A  closed  form 
analytical  solution  is  not  known  except  for  very  simple  cases,  therefore  an  approximate 
solution  based  on  numerical  methods  is  sought 

A  variety  of  methods  are  available  for  solving  the  equations  which  govern  a  turbulent 
reacting  flow.  Since  most  combustion  devices  operate  at  low  subsonic  speeds,  the 
numerical  method  used  must  be  capable  of  solving  nearly  incompressible  flows. 
Therefore,  methods  which  are  based  on  solving  density  as  a  primary  dependent  variable  are 
not  appropriate.  The  most  widely  used  methods  for  solving  such  flows  are  pressure  based 
algorithms  which  have  their  basis  in  the  SIMPLE  algorithm  of  Patankar  and  Spalding, 
1972. 

The  simulation  of  turbulent  reacting  flows  in  practical  combustion  devices  requires  use 
of  a  coordinate  system  which  is  aligned  with  the  physical  geometry  of  the  device.  The 
governing  equations  are  therefore  written  in  terms  of  general  curvilinear  coordinates 
(Thompson,  et.  al.,  1974). 

In  solving  the  flow  equations  in  curvilinear  coordinates,  one  can  solve  for  the  staggered 
Cartesian  velocity  components  as  was  done  in  the  original  SIMPLE  algorithm.  Various 
researchers  have  used  this  method  with  some  success  (Shyy,  Tong,  and  Correa,  1985; 
Correa  and  Shyy,  1987;  Shyy  and  Braaten,  1988).  It  has  recently  been  learned  however, 
that  the  convergence  rate  of  such  methods  diminishes  as  the  grid  lines  increasingly  deviate 
from  a  Cartesian  like  coordinate  system  (Vanka,  et.  al.,  1989).  Thus,  the  solution  of  the 
flow  in  such  devices  as  reverse  flow  annular  combustors  would  not  be  well  convergent. 
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Other  methods  applicable  for  solving  the  flow  in  complex  geometries  have  been 
developed.  Maliska  and  Raithby,  1984  proposed  a  method  whereby  both  Cartesian 
velocity  components  are  stored  on  each  face.  The  convergence  problems  for  reverse  flow 
geometries  are  eliminated  but  more  variables  must  be  stored  and  more  equations  must  be 
solved. 

Several  researchers  have  proposed  methods  which  solve  grid  oriented  velocities  as  the 
dependent  variables  of  the  momentum  equations.  Karki  and  Patankar,  1988  proposed  a 
method  for  solving  covariant  velocity  components  (velocities  which  are  parallel  to  the 
general  curvilinear  coordinates).  The  algorithm  involved  a  decoupled  solution  procedure  in 
which  the  continuity  equation  was  replaced  with  a  pressure  correction  equation.  The 
SIMPLER  algorithm  (Patankar,  1980)  was  used  to  correct  pressure  and  a  SIMPLE  like 
correction  equation  was  used  to  correct  the  covariant  velocities.  Joshi  and  Vanka,  1990 
proposed  a  method  based  on  contravariant  velocities  as  the  dependent  variables  (velocities 
which  are  perpendicular  to  the  coordinate  lines).  The  method  did  not  employ  a  pressure 
correction  equation,  but  rather  retained  the  continuity  equation  in  its  primitive  variable 
form.  A  coupled  solution  of  the  momentum  and  continuity  equations  (Vanka,  1986)  was 
used  to  solve  for  contravariant  velocities  and  pressure.  Demirdzic,  et.  al.,  1990  also 
proposed  a  method  based  on  contravariant  velocities  as  dependent  variables.  The  method 
involved  a  decoupled  solution  procedure  and  employed  a  pressure  correction  equation  to 
correct  for  the  contravariant  velocities  and  pressure.  One  drawback  to  procedures  which 
are  based  on  solving  for  grid  oriented  velocities  is  that  they  require  the  calculation  of  grid 
sensitive  curvature  terms  (Peric,  et.  al.,  1988).  Also,  since  the  methods  employ  a 
staggered  grid,  more  than  one  set  of  control  volumes  must  be  used  in  the  discretization  of 
the  governing  equations. 

Recent  research  in  solving  for  the  flow  in  complex  geometries  has  led  many  researchers 
to  consider  non-staggered  or  collocated  grids  for  velocities  and  pressure.  Rhie  and  Chow, 
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1983  developed  a  method  appropriate  for  solving  the  flow  equations  on  a  non-staggered 
grid.  A  method  which  avoided  the  well  known  checkerboard  split  in  the  pressure  field. 
This  method  has  since  been  extended  to  solve  for  much  more  complex  flows  (Rhie  and 
Stowers,  1988;  Rhie  and  Stowers,  1989;  Rhie  and  Syed,  1990).  Further  work  by 
Majumdar,  1988  and  Acharya  and  Moukalled,  1989  have  eliminated  some  early  problems 
with  the  original  method  proposed  by  Rhie.  These  problems  involved  relaxation 
dependencies  and  failure  to  satisfy  integral  mass  conservation.  Further  research  in  this  area 
has  been  conducted  by  Miller  and  Schmidt,  1988,  Thiart,  1990,  and  Kobayashi  and 
Pereira,  1991. 

In  the  present  work,  a  computer  program  appropriate  for  the  numerical  simulation  of 
turbulent  reacting  flow  has  been  developed.  The  simulation  is  based  on  the  solution  of 
mass  averaged  equations  on  a  non-staggered  grid.  Suitable  models  for  turbulence  and 
combustion  have  been  incorporated  into  the  numerical  method  proposed  by  Rodi, 
Majumdar,  and  Schonung,  1987. 

In  the  following  report,  the  governing  equations  for  a  turbulent  reacting  flow  are  first 
treated.  The  presentation  of  the  k-e  turbulence  model  is  followed  by  a  description  of 

models  for  diffusion  and  premixed  flames.  The  numerical  method  is  next  described  and  is 
followed  by  an  evaluation  in  several  different  flows.  Finally,  a  summary  of  the  work  done 
is  given  along  with  conclusions  and  recommendations  for  the  use  of  the  presently  proposed 
solution  algorithm. 
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CHAPTER  2:  GOVERNING  EQUATIONS 


The  fundamental  equations  which  govern  the  distribution  of  flow  describing  variables 
are  based  on  the  conservation  laws  for  mass,  momentum,  and  energy.  The  present  model 
of  turbulent  reacting  flow  seeks  to  obtain  a  solution  of  the  time  mean  steady  form  of  these 
equations. 

2.1  Laminar  Flow  Equations 


The  flow  of  a  laminar,  Newtonian  fluid  is  governed  by  the  familiar  Navier-Stokes 
equations.  In  Cartesian  x-y  coordinates,  these  equations  are  given  as  follows: 


S<P»)+^(Pv)  =  0  (2.1) 

CP«u)  +  (P«v)  =  -  5^+ (2-2) 

J;(P“v)+J^(pw)  =  -|E+^  +  ^  (2.3) 

<2-4> 

%  =  M2!-s]  <2-5> 

fdu  3vl 

^y-V'PL5y  sj  (2-6) 


|r(puho) +|t(pvi>o) =|rtea^f] +|;tea|f]  •  h°  <2-  7> 

i=l 
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*0  =  / 


CpdT  +  | 


(2.  8) 


d_ 

dx 


<PuY i>  +  ly  <*YI>  4  [fe  + 1  [fe  w\  +  <2'  9> 


8yLSc  dy 


The  present  model  is  also  applicable  for  axisymmetric  flows.  The  governing  equations  in 
r-z  coordinates  are  given  in  appendix  A.  A  detailed  derivation  of  the  above  equations  may 
be  found  in  Burmeister,  1983  or  in  Kuo,  1986. 

The  present  model  assumes  that  Fourier's  law  of  heat  conduction  is  valid  and  that 
specie  concentrations  are  sufficiently  dilute  such  that  Fick's  law  of  mass  diffusion  may  be 
used.  The  transport  properties  which  these  laws  introduce  into  the  model  are  taken  to  be 
isotropic.  There  is  thus  no  preferred  direction  for  diffusive  processes. 

It  is  assumed  that  Soret  and  Dufour  effects  are  not  present  in  the  flowfield  (Libby  and 
Williams,  1980).  Also,  viscous  heating  and  radiation  effects  are  assumed  to  be  negligible. 
The  former  assumption  is  valid  provided  that  the  flow  Mach  number  is  not  large. 

The  solution  of  the  above  equations  is  valid  only  for  laminar  flows.  Due  to  the  random 
and  highly  fluctuating  nature  of  a  turbulent  flow,  a  time  dependent  spatially  resolved 
solution  of  the  flowfield  (direct  numerical  simulation)  becomes  difficult  to  obtain  since  the 
length  and  time  scales  of  the  turbulence  are  extremely  small.  Such  a  direct  numerical 
simulation  would  require  an  exorbitant  amount  of  effort  for  even  geometrically  simple 
flows  (Moin  and  Kim,  1982).  A  stochastic  description  of  the  flow  is  thus  sought.  Rather 
than  seek  a  solution  of  the  time  dependent  properties  describing  a  turbulent  flow,  the 
averaged  flow  variables  are  the  desired  result  in  a  stochastic  description  of  a  turbulent  flow. 
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2.2  Turbulent  Flow  Equations 


The  equations  which  govern  the  distribution  of  flow  properties  in  a  turbulent  flow  are 
obtained  by  an  averaging  process.  Both  Reynolds  averaging  and  Favre  averaging  have 
been  used  in  the  description  of  turbulent  flows.  The  present  model  of  turbulent  reacting 
flow  employs  Favre  averaging. 

2.2.1  Favre  Averaging 


Favre  averaging  or  mass  weighted  averaging  is  recommended  for  variable  density 
flows  (Bilger,  1975).  Such  an  averaging  process  involves  decomposing  the  flow  variables 
into  a  suitable  mean  and  fluctuating  component: 


where, 


([>  =  <{>+<(>" 


(2. 10) 


(2.11) 


and  <)>  is  any  flow  variable  except  pressure,  viscosity,  or  density.  In  using  mass  weighted 
averaging,  the  pressure  and  density  are  decomposed  into  time  mean  and  fluctuating 
quantities,  e.g.. 


where, 


P=  P  +P 


(2. 12) 
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p  = 


(2. 13) 


The  fluctuations  in  molecular  viscosity  are  neglected.  After  decomposing  each  flow 
variable  into  its  appropriate  mean  and  fluctuating  component,  the  governing  equations  are 
time  averaged.  The  equations  which  govern  the  distribution  of  the  mean  flow  variables  are 
thus  given  as  follows  (Anderson,  et.al.,  1984): 


Jr  (pa®  +  ^  (PTO)  =  - -  pu"u"  >+  ^  -  pu  v ) 

S<p!TC)+^<PW)  =  -i |  +  J^(rxy-  pu pVV) 


dx  dy 
d  v"  d  u" 


2 

~du 

3  v 

2 

^XX  —  3  P 

L2ar 

dy- 

+  3  ^ 

2 

03v 

3  u 

2 

Tyy  -  3  JX 

LV 

3x_ 

+  3P 

dy  ’  dx 

IF  d^r 


-  [du  dv"|  ^ u  " 

Txy- V^Ldy  +  dxJ  +  M{  dy 


+  +  hQ) 


+  >C)]-|(P.v)-|(P*v)-y  *.h? 

I^i 


(2.  14) 

(2. 15) 
(2.  16) 

(2. 17) 

(2. 18) 

(2.  19) 

(2.  20) 


Jr(pori)+|?(P^i)=|  v 


(2.21) 
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4 

dy 


Sc  S?*1 


+  X  ) 


^pu-YD-^pv-Yi'-H  6* 


The  use  of  mass  weighted  averaging  instead  of  Reynolds  averaging  for  variable  density 
flows  results  in  fewer  fluctuation  terms  (e.g.  -  pu"u" )  and  makes  the  continuity  equation 

exact  (i.e.  no  terms  involving  fluctuating  quantities).  Use  of  mass  weighted  averaging 
does  however  give  rise  to  some  additional  terms  that  would  not  be  present  for  Reynolds 
averaging.  The  mean  of  the  fluctuations  such  as  in  equation  2.17  are  present  because  the 
fluctuation  is  not  one  in  time  only,  but  is  rather  a  mass  weighted  fluctuation: 

—  -  P’4>' 

<>"  = -  (2.22) 

P 

These  terms  are  small  based  on  an  order  of  magnitude  analysis  (Anderson,  et.al.,  1984) 
and  are  neglected  in  the  present  model. 

An  important  consequence  of  the  use  of  a  statistical  description  of  the  flow  is  that  the 
averaged  equations  are  not  closed.  There  are  more  unknowns  than  equations.  The 
averaged  equations  govern  the  distribution  of  the  mean  flow  variables  ( u,  v,  P,  p,  etc.  ). 

No  equations  are  available  which  govern  the  distribution  of  the  mean  fluctuations  (e.  g. 
pu"u"  ,  pu"v"  ,  etc.).  Such  a  dilemna  is  known  as  the  turbulence  closure  problem. 

A  description  of  the  mean  fluctuations  in  terms  of  the  mean  flow  variables  is  known  as 
a  turbulence  model.  Many  turbulence  models  are  available  for  the  closure  of  the  governing 
equations.  All  of  the  turbulence  models  developed  to  date  lack  universal  generality.  A 
turbulence  model  which  exhibits  generality  would  be  equally  applicable  for  a  broad  range 
of  flows  (e.g.,  a  free  shear  flow  as  well  as  a  confined  recirculating  flow).  The  available 
turbulence  models  differ  in  their  complexity  as  well  as  in  their  generality.  The  simplest  but 
more  restrictive  turbulence  models  are  the  so  called  zero  equation  models  (Anderson,  et.al.. 
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1984).  Prandtl's  mixing  length  model  is  one  of  the  most  successful  of  the  zero  equation 
models.  At  the  other  end  of  the  spectrum  are  the  Reynolds  stress  models  which  solve 
partial  differential  equations  for  the  mean  fluctuations  (also  called  Reynolds  stresses) 
(Speziale,  1991). 

The  turbulence  model  used  in  the  present  study  is  the  standard  k-e  turbulence  model. 
The  k-e  model  is  based  on  the  concept  of  an  eddy  viscosity.  The  governing  equations  are 
closed  by  describing  the  mean  fluctuations  in  terms  of  the  gradients  in  the  mean  flow  field. 

2.2.2.  Boussinesq  Approximation 

The  Boussinesq  approximation  (Anderson,  et.  al.,  1984)  is  invoked  to  relate  the 
Reynolds  stresses  to  the  gradients  of  the  mean  flow  field: 


Similar  expressions  are  used  for  the  other  turbulent  fluxes  (e.g.  pu'Tpy  ,  etc.).  Use  of 

the  Boussinesq  approximation  is  analogous  to  the  manner  in  which  the  viscous  stresses  in  a 
laminar  flow  were  related  to  the  velocity  gradients.  The  quantity  pT  is  referred  to  as  the 
turbulent  or  eddy  viscosity,  it  is  determined  in  accordance  with  the  standard  k-e  turbulence 
model, 


pT  = 


C„pk2 

e 


(2.24) 
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where  k  is  the  turbulent  kinetic  energy  and  e  the  turbulent  dissipation  rate.  The  details  of 
the  turbulence  model  are  treated  in  the  next  chapter. 

With  the  incorporation  of  the  Boussinesq  approximation  and  the  determination  of  the 
eddy  viscosity  through  knowledge  of  k  and  e,  the  equations  are  thus  closed. 


Jx  (pu)  +  (pv)  =  0 


(2. 25) 


^  (P®»)  +  ^  (Pvu>  =  -  H  + [reff  |i]  +  [reff 

^  (PUV)  +  ^  (PVV)  =  - +  J;  [reff  ^  [reff 


reg  dh0 

Pry  9x 


(Oihi 


(2. 26) 

(2. 27) 


(2.  28) 


Jx  (puYi)  +  (pv^i)  = 


Sc?  dx 


+  £ 
dy 


IkaJ[i 

Scj  dy. 


+  ODj 


(2. 29) 


The  turblent  flow  equations  for  axisymmetric  geometries  are  given  in  appendix  B.  The 
above  equations  are  used  in  the  present  study  to  model  turbulent  reacting  flows.  It  has  been 
assumed  that  the  cross  derivative  terms  in  the  momentum  equations  (  e.g.  [reff 

etc. )  are  negligible.  Such  terms  are  not  present  for  uniform  viscosity,  incompressible  flow 
and  are  here  assumed  to  be  small.  Also,  the  turbulent  Prandtl  and  Schmidt  numbers  are 
assumed  constant.  Recent  work  by  Sturgess  et.  al.,  1984  has  revealed  that  better  results 
may  be  obtained  by  allowing  these  quantities  to  vary.  However,  the  functional  relations  for 
these  quantities  are  not  certain  and  are  seemingly  problem  dependent. 
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2.3.  General  Conservation  Principle 


All  of  the  governing  equations  presented  follow  a  general  conservation  principle.  Three 
physical  processes  are  responsible  for  distributing  the  flow  variables  throughout  the 
domain  of  interest:  convection,  diffusion,  and  generation/depletion  through  a  source  or 
sink.  The  general  conservation  equation  which  accounts  for  these  processes  is  as  follows: 


reff 


sur, 

axJ+ayL1 


reff^  +  S<t>(x,y)  (2.30) 


T  +  — 7 
Prf  Pr£ 
L  T 


(2.  31) 


where  PrJ*  and  Pr£  are  laminar  and  turbulent  Prandtl  numbers  of  the  variable  <().  The 


governing  equations  for  the  various  flow  variables  are  thus  seen  to  differ  only  in  the  values 
of  diffusive  exchange  coefficient  (reff)  and  source  term  (S^). 


2.4.  Coordinate  Transformation 


The  equations  which  govern  the  distribution  of  the  flow  variables  are  transformed  to  a 
general  curvilinear  coordinate  system  so  as  to  facilitate  study  of  flows  in  complex 
geometries.  The  coordinate  transformation  involves  changing  the  independent  variables  of 
the  governing  equations.  Rather  than  describing  the  flow  geometry  in  x-y  coordinates,  it  is 
described  in  general  coordinates  (£,1)).  Figure  1  is  an  example  of  how  the  general 
coordinate  system  (£,,t\)  would  appear  in  the  Cartesian  coordinate  system  (x,y). 

The  geometries  of  practical  combustion  systems  are  such  that  they  are  not  readily 
described  by  a  Cartesian  coordinate  system.  To  avoid  use  of  a  "  stair-stepping "  procedure 
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(Rhode,  et.  al.,  1982),  a  coordinate  transformation  is  invoked  such  that  the  coordinates  of 
the  computational  grid  are  coincident  with  the  physical  boundaries  of  the  given  geometry. 
Use  of  a  coordinate  transformation  greatly  enhances  the  accuracy  of  the  boundary  condition 
specification  in  comparison  with  the  alternative,  i.e.  stair-stepping.  The  latter  procedure 
requires  the  boundary  conditions  to  be  interpolated  to  appropriate  grid  points.  If  the 
gradients  of  the  dependent  variables  are  large  in  the  vicinity  of  a  boundary  which  is  defined 
via  the  stair-stepping  procedure,  the  boundary  conditions  are  particularly  inaccurate 
(Thompson,  et.  al.,  1974). 

The  transformation  of  coordinates  is  implemented  by  considering  that  the  Cartesian 
coordinates  (x,y)  are  functions  of  the  general  curvilinear  coordinates  (^,r|),  whence, 

x  =  fj(4,T|) 

y  = 

The  transformation  must  be  one  to  one.  A  given  point  (x,y)  can  correspond  to  only  one 
point  in  the  transformed  plane  (£,r|)  (Fletcher,  1988b).  All  derivatives  in  x,y  are 
subsequently  expanded  as  functions  of  ^,T|  through  use  of  the  chain  rule  of  calculus: 


3<3> _  3<f>  d<)>  1 

^  L^yi1 '  dny*\|J 


dtj)  _  5(J) 

dy  L  a^XT1 


(2.  32) 

(2.  33) 


where  J  =  x^  -  x^. 

The  transformation  enforced  in  the  present  method  is  such  that  the  size  of  a  cell  in  the 
transformed  plane  is  unity  (Al;  =  Atj  =  1).  In  general,  an  analytical  expression  is  not 
available  for  determining  the  transformation  and  so  it  is  determined  numerically 
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(Thompson,  et.  al.,  1974).  The  present  method  is  appropriate  for  use  with  algebraic  and 
differential  equation  methods  of  defining  the  transformation  (grid  generation). 

Upon  expanding  all  the  derivatives  as  functions  of  t,  and  tj  and  by  rearranging,  the 
general  conservation  equation  may  be  written  as  follows: 


|(pu,)+^(pv«=|[rq,1|]+|)[rq22|] 

+i[rqi2^]+^[r,2iS] 

+msns,Ti) 

(2.  34) 

U  =  uyT1-vx11 

(2.  35) 

V  =  vx^  -  uy^ 

(2.  36) 

qn  =  i  J  |(yT12+V) 

(2.  37) 

922  =  1  j  |  ( y^2  +  x42 ) 

(2.  38) 

Q2i = qi2 = -  ni  + ) 

(2.  39) 

1 J 1  =  ( xjov  Vs ) 

(2. 40) 

The  quantities  U  and  V  are  contravariant  velocities.  The  contravariant  velocity  U  is 
perpendicular  to  coordinate  lines  of  constant  £  and  is  actually  a  volume  flow  rate.  The 
contravariant  velocity  V  is  perpendicular  to  coordinate  lines  of  constant  r\  and  is  also  a 
volumetric  flow  rate.  For  an  incompressible  flow,  these  quantities  are  conserved.  The 
transformation  metrics  corresponding  to  a  transformation  from  r,z  to  coordinates  are 
different  than  those  given  above.  These  transformation  metrics  are  given  in  appendix  C. 

The  governing  equations  in  the  transformed  space  may  be  written  in  several  different 
forms.  The  governing  equations  have  been  written  in  conservation  law  form  in  the  present 
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implementation  since  this  form  leads  naturally  to  conservative  discretizations  (Anderson,  et 
al.,  1984). 

2.5.  Boundary  Conditions 

The  prescription  of  a  well  posed  mathematical  problem  requires  specification  of 
appropriate  boundary  conditions.  Since  the  governing  equations  are  written  as  functions  of 
the  independent  variables  ^  and  T|,  the  boundary  conditions  are  imposed  in  the  transformed 

plane.  In  the  present  method,  four  types  of  boundary  conditions  have  been  utilized:  wall, 
symmetry,  inlet,  and  outlet 

2.5.1.  Wall  Boundary  Condition 

All  solid  boundaries  are  treated  as  no  slip,  no  penetration  boundaries.  The 
contravariant  velocities  (U  and  V)  as  well  as  the  Cartesian  velocities  are  equal  to  zero  at  the 
walls. 

The  walls  are  also  treated  as  being  adiabatic.  The  derivative  of  the  stagnation  enthalpy 
in  the  direction  perpendicular  to  the  wall  is  equal  to  zero.  For  example,  if  a  wall  lies  along 
a  coordinate  of  constant  T|,  then  the  following  is  specified, 

i:=0  (2-41) 

Such  a  boundary  condition  is  only  true  if  the  grid  is  orthogonal  at  the  boundary.  If  not, 
some  error  is  introduced  at  the  boundary.  The  error  will  be  proportional  to  the  extent  of 
non-orthogonality  at  the  wall.  Thompson  gives  a  boundary  condition  implementation 
which  is  appropriate  for  boundary  non-orthogonal  grids  (Thompson,  et.  al.,  1974). 
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A  similar  treatment  is  used  for  the  fuel  fraction  at  a  wall,  thus  treating  all  solid 
boundaries  as  non-catalytic.  The  boundary  conditions  for  the  turbulence  quantities  k  and  e 
are  imposed  by  wall  functions.  Details  of  the  wall  function  implementation  are  given  in  the 
next  chapter. 

2.5.2.  Symmetry  Boundary  Condition 

A  zero  first  derivative  is  specified  perpendicular  to  any  boundary  designated  as  a 
symmetry  boundary.  The  zero  derivative  is  specified  for  the  Cartesian  velocity  components 
as  well  as  the  scalar  quantities.  The  contravariant  velocities  perpendicular  to  any  symmetry 
boundary  are  specified  as  zero.  The  contravariant  velocities  parallel  to  a  symmetry  are 
treated  in  the  same  manner  as  scalar  quantities. 

2.5.3.  Inflow  Boundary  Condition 

The  present  model  of  turbulent  reacting  flow  treats  the  inlet  as  a  Dirichlet  boundary 
condition.  Values  of  velocities  and  scalars  are  specified  at  the  inlet.  A  uniform  distribution 
of  the  flow  variables  is  taken  as  the  default  specification.  Parabolic  and  other  non-uniform 
distributions  of  flow  variables  can  also  be  used.  Currently,  the  model  does  not  account  for 
pressure  boundary  conditions,  thus  phenomena  such  as  choking  cannot  be  treated.  The 
inlet  specification  is  best  suited  for  incompressible  and  subsonic  compressible  flows. 
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2.5.4.  Outflow  Boundary  Condition 

A  fully  developed  flow  condition  is  specified  for  all  variables  at  any  boundary  that  is 
designated  as  an  outflow.  If  an  outflow  lies  along  a  coordinate  line  of  constant  then  the 

following  would  be  specified  at  the  outflow  plane: 


The  contravariant  velocities  parallel  to  the  outflow  plane  are  additionally  specified  as  being 
equal  to  zero.  The  outflow  boundary  condition  has  no  affect  on  the  solution  domain  if  the 
outflow  plane  is  properly  specified.  With  the  given  discretization  scheme  (see  Ch.  4),  all 
information  at  the  outlet  plane  is  cut  off  from  the  rest  of  the  domain  when  the  cell  Peclet 
number  is  greater  than  two  (Patankar,  1980).  Thus,  the  outflow  plane  is  properly  specified 
if  the  flow  in  that  region  is  expected  to  be  locally  "one-way"  (i.e.  no  recirculation)  and  if 
the  convective  effects  are  dominant  (Peceu  >  2). 

2.6.  Thermodynamic  and  Transport  Properties 

The  prescription  of  thermodynamic  and  transport  properties  is  required  to  close  the 
system  of  equations  governing  the  flow.  The  specific  heat,  molecular  weight,  molecular 
viscosity,  specific  heat  ratio,  and  density  all  need  to  be  related  to  the  primary  dependent 
variables  (u,  v,  P,  h0,  Ypu,  etc.) 

The  fluid  density  is  determined  via  the  equation  of  state.  It  is  assumed  that  the  fluid 
obeys  the  ideal  gas  law,  hence: 
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(2. 43) 


_  PW 
p  "  RUT 

where  W  is  the  mixture  molecular  weight  and  Ru  is  the  universal  gas  constant 

The  molecular  weight  is  a  function  of  the  specie  mass  fractions.  If  Dalton's  model  of 
partial  pressures  is  taken  to  describe  the  reacting  mixture,  then  the  mixture  molecular 
weight  is  given  as  follows: 


(2.44) 


The  constant  pressure  specific  heat  is  determined  by  the  specie  mass  fractions  and  by 
the  temperature.  Following  the  treatment  used  in  the  NASA  thermochemistry  code 
(Gordon  and  McBride,  1976),  the  specific  heat  for  a  particular  specie  is  given  by  the 
following: 


Cpi  —  ( ai  +  <*2^  +  ^T^ )  (2. 45) 

The  mixture  specific  heat  is  defined  as  a  mass  weighted  average  of  the  specie  mass 
fractions: 


N 

Cpmix  =  ^  YjCpj 

i=l 


(2. 46) 


The  specific  heat  ratio  is  taken  to  be  a  constant  in  the  present  model.  The  specific  heat  ratio 
for  air  (y=1.4)  is  assumed  to  be  an  accurate  representation  of  the  actual  value  of  the  gas 
mixture. 
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In  the  present  model  for  turbulent  reacting  flow,  the  molecular  viscosity  is  assumed  to 
be  a  function  of  only  temperature.  The  molecular  viscosity  is  determined  from  the 
Sutherland's  law  expression  for  air  (White,  1974): 


(X  rT]3/2T0+S 

mThoJ  t  +  s 


(2.47) 


A  more  rigorous  treatment  of  the  molecular  viscosity  would  treat  the  various  species 
present  in  the  flow.  Such  a  detailed  treatment  of  the  molecular  viscosity  is  not  considered 
necessary  since  the  turbulent  "eddy"  viscosity  is  much  larger  (Bilger,  1980).  The  exchange 
coefficient  appearing  in  the  governing  equations  is  thus  primarily  determined  by  the 
turbulent  viscosity.  For  the  calculation  of  laminar  reacting  flows,  a  more  detailed  treatment 
of  the  molecular  viscosity  would  be  beneficial. 
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CHAPTER  3:  TURBULENCE  AND  COMBUSTION  MODELS 


"  Mathematical  models  are  tools  used  in  engineering  and  science  to  predict  functional 
relationships  between  certain  input  and  output  functionals."  (Magnussen  and  Hjertager, 
1976).  In  simulating  turbulent  reacting  flows,  mathematical  models  are  employed  in  order 

to  close  the  system  of  equations  which  governs  the  mean  flowfield.  The  desired  "  output " 
of  these  models  are  Reynolds  stresses  ( -  pu"v"  ,  etc.),  turbulent  fluxes  of  scalar  quantities 

(e.g.  pu"  Y  ),  and  mean  chemical  reaction  rates  (  d) ).  The  input  variables  of  the  models 

are  properties  of  the  mean  flowfield  (i.e.  u,  v,  YFU,  etc.  ).  The  mean  properties  are 
governed  by  the  fundamental  conservation  laws.  The  turbulence  and  combustion  models 
thus  relate  the  turbulent  fluctuation  quantities  to  the  primary  dependent  variables. 

The  present  simulation  of  turbulent  reacting  flows  uses  the  standard  k-e  turbulence 
model  with  wall  functions  to  facilitate  closure  of  the  governing  equations  (Launder  and 
Spalding,  1974).  Two  simple  models  of  combustion  are  given  in  the  present  simulation:  a 
turbulent  premixed  flame  model  and  a  diffusion  flame  model.  The  turbulent  premixed 
flame  simulation  is  based  on  an  EBU  model  proposed  by  Magnussen  and  Hjertager,  1976. 
The  diffusion  flame  model  is  based  on  the  assumption  of  infinitely  fast  kinetics  and  is  of 
the  assumed  pdf  variety. 

3.1.  Turbulence  Model 


The  objective  of  the  turbulence  model  is  to  relate  the  turbulent  fluxes  (e.g.  -  pu"v" , 

etc.)  to  the  mean  flow  variables  (u,  v,  etc. ).  As  was  stated  in  the  previous  chapter, 
there  are  in  existence  numerous  procedures  for  obtaining  turbulence  closure.  The  various 
methods  range  from  the  simple  mixing  length  models  to  the  more  complex  Reynolds  stress 
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models.  The  available  turbulence  models  all  differ  in  their  complexity  and  generality.  For 
the  present  effort,  it  was  desired  to  use  a  computationally  efficient  turbulence  model.  A 
model  which  exhibited  reasonable  generality  without  requiring  an  exhaustive  computational 
effort.  The  standard  two  equation  k-e  turbulence  model  with  wall  functions  is  used  in  the 
present  simulation. 

The  standard  k-e  model  is  based  on  the  concept  of  a  turbulent  or  "eddy"  viscosity.  The 
Reynolds  stresses  are  written  in  terms  of  the  mean  flow  gradients  by  invoking  the 
Boussinesq  assumption: 


ti  tt 

-  pUjUj  =pT 


2  ~  T  3uk 
'  3  5y  l_^T  8xT 


(3. 1) 


A  similar  relation  is  used  for  expression  of  other  turbulent  fluxes.  The  use  of  the 
Boussinesq  approximation  shifts  the  burden  of  closure  from  one  of  describing  -  pu"v"  ,  to 

one  of  determining  the  eddy  viscosity  jij.  The  modelling  of  the  eddy  viscosity  is  achieved 

by  analogy  with  the  kinetic  theory  of  gases. 

The  expression  for  the  viscosity  of  a  gas  as  obtained  from  kinetic  theory  is  as  follows 
(Burmeister,  1983): 


vAmn 

n=— 


(3.2) 


where, 

v  =  average  particle  velocity 
A,  =  mean  free  path  of  particles 
m  =  particle  mass 

n  =  number  of  particles  per  unit  volume 
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In  an  analogous  manner,  the  eddy  viscosity  is  taken  to  be  directly  proportional  to  turbulent 
length  and  velocity  scales, 


M-t  ~  pV’pLj 


(3.  3) 


The  k-e  model  uses  the  square  root  of  turbulent  kinetic  energy  (k1/2)  as  the  characteristic 
velocity.  The  turbulent  length  scale,  Lj  is  determined  as  a  function  of  the  turbulent 
dissipation  rate,  e  (Launder  and  Spalding,  1974): 


L  j  = 


CDk3/2 

e 


(3.4) 


The  resulting  expression  for  the  eddy  viscosity  is  therefore: 


(3.  5) 


The  k-e  turbulence  model  determines  both  the  turbulent  velocity  and  length  scales  from 
transport  equations: 


a$L  K  dnJ  dnl 

+  1  JIPk-pel  Jl 


(3.6) 


|(pUe)  +  JW| 


'  de‘ 

d 

del 

r- qi1  iil 

+ - 

dn 

E  ,22 

23 


(3.7) 
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L  driJ  dq 
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^^t{2[|]2  +  2[|]2  +  [|+^}  (3.8) 


It  is  evident  from  the  expression  for  eddy  viscosity  that  it  is  a  property  of  the  flow  and 

not  the  fluid.  Different  flow  geometries  and  inlet  conditions  produce  varied  distributions  of 
k  and  e  and  thus  jit. 

The  k-e  model  is  a  reasonably  simple  means  of  closure  as  regards  computational  effort. 
Only  two  additional  differential  equations  are  introduced  into  the  solution  procedure.  The 
price  paid  for  the  simplicity  is  lack  of  universality.  The  k-e  model  treats  the  turbulence  as 
having  only  a  single  characteristic  velocity  and  length  scale  respectively.  Such  a 
description  is  not  valid  for  all  flows  (Burmeister,  1983).  Two  additional  constraining 
assumptions  incorporated  into  the  k-e  turbulence  model  are  that  the  model  is  based  on  local 

isotropy  and  the  presence  of  fully  developed  turbulence.  The  latter  assumption  requires 
near  wall  regions  to  be  treated  in  a  special  manner  since  the  flow  in  the  immediate  vicinity 
of  a  solid  boundary  is  not  fully  turbulent  (Schlichting,  1979).  The  near  wall  region  may  be 
treated  by  employing  an  alternate  form  of  the  k-e  model  known  as  a  low  Reynolds  number 
k-e  model  (Patel,  et.  al.,  1985).  The  low  Reynolds  number  k-e  models  employ  suitable 
modifications  of  the  transport  equations  for  k  and  e  in  near  wall  regions.  In  using  these 
models  it  is  necessary  to  resolve  the  near  wall  region  with  the  chosen  computational  grid. 
An  alternative  procedure  is  to  use  wall  functions  as  a  means  of  linking  the  near  wall  flow  to 
the  fully  turbulent  core  (Launder  and  Spalding,  1974). 
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Wall  functions  assume  that  the  flow  near  a  solid  boundary  is  locally  a  Couette  flow 
(Burmeister,  1983).  With  this  assumption,  the  component  of  velocity  parallel  to  the  wall  is 
described  by  a  logarithmic  profile  (Schlichting,  1979): 

u+  =  ~  In  (Ey+)  (3.9) 

Jv 

K  =  0.4 
E  =  9.025 


It  is  further  assumed  that  the  production  and  dissipation  of  turbulent  kinetic  energy  are 
equal  in  the  near  wall  region  (Hwang  and  Liou,  1991),  thus  implying  the  following: 


k  = 


C 


P 


(3. 10) 


The  logarithmic  velocity  profile  is  thus  expressed  as  follows: 


_  C^kl/2, 

u  =— 14 - In 

K 


EyCg^klgp' 

P 


(3. 11) 


A  straightforward  manner  of  implementing  the  wall  functions  for  velocity  is  to  replace  the 
near  wall  velocity  with  that  obtained  from  the  log-law  profile.  (This  is  done 
computationally  through  suitable  source  term  linearization  (Patankar,  1980)  ).  In  the 
context  of  a  geometry  described  by  general  curvilinear  coordinates,  this  can  be  a  tedious 
process  since  the  log-law  profile  is  utilized  for  the  component  of  velocity  which  is  parallel 
to  the  wall.  In  general  curvilinear  coordinates,  this  velocity  will  be  comprised  of  both 
Cartesian  velocity  components.  Rather  than  forcing  the  near  wall  velocity  to  be  equal  to  the 
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log-law  value  as  proposed  above,  the  near  wall  eddy  viscosity  is  suitably  modified  so  as  to 
reflect  the  resulting  turbulent  shear  stress, 


pKC,1^kV2yn 

"ln  rECu1/4k1/2Pyp 

L  1 1 


C3. 12) 


where  yp  is  the  distance  from  the  wall  to  the  center  of  the  cell  nearest  the  wall.  In  using  the 

wall  function  method,  it  is  important  to  note  that  the  cell  nearest  the  wall  must  be  such  that 
yp  lies  within  the  fully  turbulent  core. 

The  near  wall  affect  on  turbulent  kinetic  energy  is  enforced  by  writing  the  production 
and  dissipation  terms  in  the  k  transport  equation  in  their  integrated  form  according  to  the 
log-law  velocity  profile.  Therefore, 


^  u [EC"1/4kl/2pyn 


(3.  13) 


yp 


e  =  JJ  e  dx  dy 

C  3/4 

=  ^U__k3/2  ln 

KyD 


EC„i%i/2pyp 

M- 


(3.  14) 


The  transport  equation  for  k  is  solved  at  the  near  wall  cell  but  the  production  and 
dissipation  terms  are  modified  as  given  above  (Launder  and  Spalding,  1974). 

The  near  wall  effect  on  the  turbulent  dissipation  rate,  e  is  implemented  in  the  present 
method  by  fixing  the  near  wall  value  of  e.  Given  the  definition  of  the  turbulent  dissipation 
rate  (Launder  and  Spalding,  1974), 
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(3. 15) 


e  =  v 


du^duVj 
3xk  3xk 


the  following  approximation  is  made  in  the  near  wall,  locally  Couette  flow,  region: 


e  =  ^ 


du 

dy 


(3.  16) 


where  y  is  the  direction  normal  to  the  wall  and  u  is  the  velocity  parallel  to  the  wall.  Since 
the  mean  velocity  in  the  near  wall  region  is  taken  to  be  governed  by  a  log-law,  the  value  of 
the  turbulent  dissipation  rate  at  the  point  nearest  the  wall  is  given  as  follows: 


eP- 


Kyp 


(3. 17) 


3.2.  Premixed  Flame  Model 

The  purpose  of  the  turbulent  premixed  flame  model  is  to  provide  a  functional  relation 
for  the  determination  of  the  mean  reaction  rates  of  all  relevant  species.  A  complete 
description  of  hydrocarbon  combustion  can  require  the  consideration  of  as  many  as  83 
chemical  reactions  with  27  participating  chemical  species  (Jachimowski,  1984).  Therefore, 
27  expressions  for  d)  j  would  be  required,  one  for  each  specie. 

In  order  to  keep  the  overall  flowfield  simulation  computationally  tractable,  a  simple 
one-step  description  of  the  combustion  process  is  considered.  The  one-step  mechanism  of 
C3H8  chemistry  proposed  by  Westbrook  and  Dryer,  1981  is  used  in  the  present  simulation: 

C3H8  +  5(  02  +  3.76  N2 ) - >  3  C02  +  4  H20  +  1 8.8  N2  (3. 1 8) 
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Reaction  rate  expressions  for  the  other  species  (O2,  CO2,  etc.)  are  not  needed,  since  the 
mass  fractions  of  species  other  than  fuel  are  determined  from  the  law  of  conservation  of 
atoms.  The  number  of  carbon  atoms  entering  the  reaction  is  the  same  as  the  number  of 
carbon  atoms  after  reaction  (a  statement  of  mass  conservation),  therefore  an  expression  for 
the  mass  fraction  of  CO2  can  be  determined: 

3WC0 9 

Yco2=^72(f-YFu)  (3-21) 

Similar  expressions  exist  for  H2O  and  O2: 

4Wh?o 

Yh20— w^<f-Ypu>  (3.22) 

Yq2  =  0.2331  (  1  - f )  -  (f-Ypu)-^2  (3.23) 


The  mass  fraction  of  nitrogen  is  calculated  from  the  sum  of  all  mass  fractions  being  equal  to 
unity,  therefore: 

Yn2  =  1  •  Yrj  -  y02  -  Yh2o  -  YC02  (3.  24) 
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In  describing  the  mean  reaction  rate  (  co  )  it  is  not  sufficient  to  simply  substitute  mean 

values  of  flow  variables  into  the  reaction  rate  expression.  Such  a  procedure  neglects 
turbulent  fluctuations  of  specie  concentrations  and  temperature  and  will  lead  to  errors  of 
typically  an  order  of  magnitude  (Jones  and  Whitelaw,  1982). 

The  present  model  of  turbulent  premixed  flames  determines  the  mean  reaction  rate  as 
being  the  smaller  of  two  limiting  cases.  The  first  limiting  case  involves  a  kinetically  limited 
flame  for  which  the  decay  rate  of  turbulent  concentration  fluctuations  is  much  greater  than 
the  instantaneous  rate  at  which  the  reaction  proceeds  (Libby  and  Williams,  1980).  For  this 
limit,  the  mean  reaction  rate  does  equal  the  expression  given  above  (Zhou,  1989): 

®  ARR  =  '^f  P^+b^  ^fu ?02  (3. 25) 


=  B  exp  (3. 26) 

The  expression  for  the  mean  reaction  rate  constant  (fq)  is  further  simplified  in  the  present 
model: 


kf  =  B  exp 


(3.  27) 


This  expression  assumes  that  the  decay  rate  of  temperature  fluctuations  is  also  much  greater 
than  the  instantaneous  reaction  rate. 

The  other  limiting  case  involves  combustion  in  which  the  instantaneous  chemical 
reaction  occurs  so  fast  that  only  the  turbulent  mixing  limits  the  combustion  process.  The 
chemical  reaction  rate  is  much  greater  than  the  decay  rate  of  turbulent  fluctuations.  This 
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limit  corresponds  to  a  physically  controlled  turbulent  flame  where  the  chemical  kinetics  do 
not  effect  the  flame  shape  and  length.  An  expression  for  the  mean  reaction  rate  for  this 
limit  is  given  by  Magnussen  and  Hjertager,  1976: 


®  EBU  =  ■  CEBU  P  ^FUu 


(3. 28) 


Two  different  expressions  are  available  for  the  determination  of  the  mean  reaction  rate. 
In  the  simulation  of  turbulent  premixed  flames,  the  mean  reaction  rate  is  taken  as  equal  to 
the  smaller  of  the  two  limiting  cases  given.  Therefore, 

to  fu  = "  n“n  l  *  to  arr  I  »  I  to  ebu  I  )  (3.  29) 

3.3.  Diffusion  Flame  Model 

The  present  simulation  of  turbulent  reacting  flow  considers  two  models  for  diffusion 
flames.  The  first  model  is  for  a  laminar  flame.  The  second  model  is  applicable  for 
turbulent  diffusion  flames  and  is  based  on  the  fundamental  concepts  of  the  laminar  flame 
model. 


3.3.1.  Laminar  Flame  Model 

A  model  for  laminar  flames  is  not  required  as  regards  obtaining  a  closure  of  equations. 
Since  an  averaging  process  is  not  used,  no  turbulent  fluxes  or  mean  reaction  rates  are 
present  in  the  governing  equations.  However,  a  model  is  incorporated  in  the  present 
simulation  so  as  to  simplify  the  prediction  of  laminar  diffusion  flames. 


30 


The  quantities  of  interest  in  the  description  of  a  reacting  flow  are  specie  mass  fractions 
and  the  temperature.  Rather  than  solve  the  governing  differential  equations,  it  is  assumed 
that  the  chemical  reactions  occur  at  an  infinitely  fast  rate.  This  assumption  obviates  one 
from  having  to  solve  all  the  differential  equations  which  govern  the  distributions  of  mass 
fractions  and  temperature.  Instead,  the  knowledge  of  the  mixture  fraction  is  sufficient  for 
the  determination  of  mass  fractions  and  temperature. 

The  present  model  assumes  that  the  reaction  of  fuel  and  oxidant  occurs  irreversibly  in 
one  step: 


Fuel  +  Oxidizer - >  Product 


(3.  30) 


Given  this,  knowledge  of  only  three  mass  fractions  is  required:  fuel,  oxidizer,  and 
product. 

In  the  present  model  for  laminar  diffusion  flames,  it  is  assumed  that  the  reaction  occurs 
only  when  the  fuel  and  oxidizer  are  in  stoichiometric  proportions,  thus: 

1  kg  Fuel  +  (()  kg  Oxidizer - >  ( 1  +  <|> )  kg  Product  (3.  31) 


where  <])  is  the  stoichiometric  air  to  fuel  ratio  based  on  mass.  With  the  further  assumption 
that  the  diffusivities  of  fuel  and  oxidizer  are  equal,  the  differential  equations  governing  the 
distributions  of  fuel  and  oxidizer  mass  fractions  can  be  linearly  combined  to  yield  an 
equation  for  the  Schvab-Zeldovich  variable  *P  (Bilger,  1980): 
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(3.  32) 
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¥  =  Ypu- Y0x/ <t> 


(3.  33) 


The  normalization  of  the  conserved  variable  'P  yields  a  relation  for  the  mixture  fraction  f: 


f  = 


y-^OX 
^FU  -  ^OX 


(3.  34) 


where  Tqx  =  -1/<|>  and  'Ppy  =  1  if  no  diluent  is  present  in  either  stream.  The  mixture 
fraction  is  also  a  conserved  variable  and  is  thus  governed  by  the  following  equation: 


kpm+i 


'»a+Sr«*8 


(3.  35) 


The  assumption  of  infinitely  fast  chemistry  implies  that  fuel  and  oxidizer  cannot  exist  at 
the  same  point  and  time  in  the  flowfield.  Therefore  if  the  conserved  variable  VP  is 
calculated  to  be  positive  (implying  fuel  rich)  then  only  fuel  and  product  exist  at  that  point. 
Alternately,  if  *F  is  negative  (fuel  lean  condition)  then  only  oxidizer  and  product  are 
present. 


¥>0 

Yfu  =  ^ 
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ypr  -  1  -  Ypu 

(3.  36) 

¥<0 

Yfu  =  0 
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ypr  =  1  -  Yox 

(3.  37) 

¥  =  0 

Yfu  =0 
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II 
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ypr  =  i 

(3.  38) 
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Note  that  ¥  =  0  corresponds  to  the  stoichiometric  condition  and  hence  the  location  of  the 
flame  surface.  In  the  present  model,  the  differential  equation  for  'P  is  not  solved,  rather, 
the  mixture  fraction  is  solved  from  which  *P  is  determined  according  to  the  above  relations. 

It  is  assumed  in  the  present  model  that  the  Lewis  number  is  unity  (Pr  =  Sc)  and 
therefore  the  governing  equation  for  enthalpy  is  free  of  a  source  term  (Libby  and  Williams, 
1980).  It  is  also  assumed  that  the  combustion  process  occurs  adiabatically.  The  enthalpy 
can  therefore  be  determined  as  a  function  of  the  mixture  fraction: 


h  -  f  hFUEL  +  ( 1  -  f )  h0x 


(3.  39) 


where  hpi  jpT  and  Iiqx  are  the  enthalpies  of  the  fuel  and  oxidant  streams  respectively.  The 
mixture  temperature  is  determined  from  knowledge  of  the  enthalpy: 


h-YpuQ 

Cpmix 


(3. 40) 


where  Q  is  the  heat  of  combustion. 

3.3.2.  Turbulent  Flame  Model 

The  model  for  turbulent  diffusion  flames  assumes  that  the  reactions  occur  infinitely 
fast,  just  as  in  the  laminar  case.  In  the  case  of  a  turbulent  flow,  the  quantities  to  be 
calculated  by  the  combustion  model  are  averaged  quantities:  YFU,  Yox,  YPR,  T.  For  the 
purpose  of  simplicity,  a  one-step  irreversible  reaction  is  again  assumed  to  describe  the 
kinetics.  If  the  mean  quantities  are  determined  in  a  manner  identical  to  that  for  a  laminar 
flame,  significant  error  results.  To  do  so  would  be  to  neglect  the  presence  of  turbulent 
fluctuations  in  the  flow  variables.  Such  a  treatment  of  combustion  would  assume  that  the 
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averaged  fuel  and  oxidizer  mass  fractions  cannot  exist  at  the  same  point.  This  is  not 
observed  experimentally  (Bilger,  1980).  In  order  to  determine  the  mean  flow  variables,  the 
definition  of  an  average  is  used  (Hogg,  et.  al.,  1987): 


4>P(f)df 


0 


(3. 41) 


where  P(f)  is  the  probability  density  function. 

The  quantities  to  be  determined  are  mass  averaged  values  and  therefore  a  Favre  pdf 
(Bilger,  1980)  defines  the  mean  quantities: 


YFu=/ 


YFU=  Ypu(f)  P(f:x)  df 


Yqx  -  / 


YOX“  Yox(f)P(f:x)df 


-/ 


YPr=  YPR(f)P(f:x)df 


h  =  / 


h(f)  P(f:x)  df 


(3.42) 


(3.  43) 


(3.44) 


(3.45) 


where  P(f;x)  is  the  Favre  pdf.  The  probability  density  function  used  in  the  present 
simulation  is  a  double  delta  function  pdf  (Jones  and  Whitelaw,  1982). 
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P(f:x)  =  a  8(f -f+)  +  ( 1  -  a )  8(f  -  f ) 


(3. 46) 


a  = 


f  -  f- 

f+  .  f- 


(3.47) 


Other  researchers  have  used  a  clipped  Gaussian  pdf  (Lockwood  and  Naguib,  1975),  a  beta 
function  pdf  (Lentini,  1990),  and  a  10  points  pdf  (Chen,  et.  al.,  1991). 

The  probability  density  function  is  spatially  dependent  and  methods  have  been  used 
which  solve  for  the  pdf  prior  to  determining  the  means  of  the  other  flow  variables  (Pope, 
1990).  The  model  used  in  the  current  simulation  of  turbulent  reacting  flow  is  an  assumed 
shape  pdf  model  (Kuo,  1986). 

The  double  delta  function  is  specified  by  the  fluctuation  limits  f+  and  f" : 

f+  =  f  +  Vg  (3.48) 


f  =  ?  -  Vi 


(3.49) 


The  concentration  fluctuation  g  (where  g  =  f '2  ),  is  determined  from  the  modelled  equation 
given  by  Spalding,  1971: 


(3.  51) 
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With  the  use  of  a  double  delta  function  pdf,  a  numerical  integration  is  not  required  in 
determining  mean  quantities.  The  determination  of  mean  flow  variables  which  describe  a 
reacting  flow  are  calculated  based  on  a  knowledge  of  a: 


^fu  =  «YJu  +  (1-o)Y^j 

(3.52) 

<?ox  =  <x  +  n.a)Yox 

(3. 53) 

?PR  =  aYjR  +  (l-a)Y;R 

(3.  54) 

h  =  ah+  +  (l-a)h' 

(3.  55) 

where  Ypy,  YqX,  etc.  are  the  quantities  which  correspond  to  the  mixture  fraction  at  the 
extreme  upper  limit  of  fluctuation  (f+).  The  values  YpUS  Y^x,  etc.  are  the  property 

fluctuations  corresponding  to  the  lower  limit  of  the  mixture  fraction  fluctuation  (f*).  The 
limit  values  are  determined  in  a  manner  identical  to  that  for  a  laminar  flame  except  they  are 
based  on  a  mixture  fraction  fluctuation  (f+  or  f-).  The  mean  temperature  is  determined  from 
the  mean  enthalpy  as  follows: 


~  _  h  -  Yfu  Q 


^pmix 


(3. 56) 


In  using  a  pdf  method,  the  effects  of  concentration  fluctuation  are  accounted  for  to 
some  extent.  The  ramification  of  pdf  shape  on  the  results  has  been  studied  by  various 
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researchers  (e.g.  Chen,  et.  al.,  1991).  They  concluded  that  the  shape  of  the  pdf  does  not 
greatly  affect  the  flowfield  except  when  a  double  delta  function  pdf  is  used.  The  use  of  a 
more  descriptive  yet  computationally  tractable  pdf  is  thus  desired. 
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CHAPTER  4:  SOLUTION  PROCEDURE 


To  facilitate  the  study  of  phenomena  relevant  to  turbulent  reacting  flows,  the  governing 
conservation  equations  must  be  solved.  These  equations  are  highly  nonlinear  and  are 
strongly  coupled,  thus  precluding  the  determination  of  a  closed  form  analytic  solution. 
Rather  than  solve  for  continuous  functions  of  the  flow  variables,  the  equations  are 
discretized  and  values  of  the  flow  variables  at  only  certain  points  are  sought.  The 
discretization  process  simplifies  the  problem  from  one  of  solving  partial  differential 
equations  to  one  of  solving  algebraic  equations.  The  nonlinearity  of  the  algebraic  equations 
as  well  as  the  inter-equation  coupling  is  treated  by  the  process  of  iteration.  In  general,  the 
equations  are  solved  in  a  sequential  fashion.  While  one  flow  variable  is  solved,  the  others 
are  held  constant.  By  solving  each  equation  in  turn,  one  iteration  is  completed.  Numerous 
iterations  are  required  before  a  converged  solution  is  obtained.  Convergence  is  obtained 
when  each  flow  variable  satisfies  its  discretized  equation  and  when  the  flow  variables  no 
longer  change  from  one  iteration  to  the  next. 

A  solution  procedure  employing  this  sequential  technique  is  referred  to  as  a  segregated 
or  decoupled  solution  procedure.  Coupled  solution  procedures  have  also  been  developed 
(Vanka,  1986;  Galpin  et.al.,  1985).  The  coupling  of  all  equations  is  usually  not 
undertaken  however.  Coupled  solution  procedures  normally  couple  only  the  momentum 
and  continuity  equations,  with  the  various  scalar  equations  being  solved  in  a  sequential 
decoupled  framework. 

The  present  solution  procedure  employs  a  decoupled  solution  technique.  The 
momentum  equations  are  solved  first,  followed  by  a  pressure  correction  to  satisfy  the 
continuity  constraint.  The  scalar  equations  for  turbulent  kinetic  energy,  turbulent 
dissipation  rate,  etc.  are  solved  next,  with  an  update  of  thermodynamic  properties 
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(p,T,etc.)  and  transport  coefficients  (r)  completing  one  iteration.  In  the  discussion  which 
follows,  the  discretization  of  the  momentum  equations  will  be  treated.  The  use  of  a 
collocated  grid  for  velocities  and  pressure  and  its  ramification  on  the  discretized  momentum 
equations  will  be  discussed.  The  pressure  correction  procedure  employed  will  be  described 
followed  by  the  solution  of  the  scalar  equations.  Lastly,  a  synopsis  of  the  entire  solution 
procedure  will  be  presented. 

4.1.  Momentum  Equations 

The  momentum  equations  are  discretized  in  accordance  with  a  finite  volume  method. 
The  finite  volume  method  has  been  chosen  because  of  the  conservative  nature  of  the 
discretization  obtained  (Anderson,  et.al.,  1984;  Fletcher,  1988a).  In  discretizing  the 
momentum  equations  using  a  finite  volume  method,  the  computational  domain  is  first 
subdivided  into  finite  control  volumes  (see  Figure  2).  The  Gauss  divergence  theorem  is 
then  used  to  integrate  the  partial  differential  equation  over  a  given  finite  control  volume. 
The  resulting  discrete  forms  of  the  u  and  v  momentum  equations  are  given  as  follows: 


[  (pUu)i+1/2 j  -  (pUu)j.1/2j  ]  Ar|  +  [  (pVu)iJ+1/2  -  (pVu)ij.1/2  ]  A£,  = 

[  (Tqil^  )i+l/2j  -  (Tqnj|  )i-i/2j  ]  ATI  +  (4.  1) 

[  )ij+l/2  -  >i  j-1/2  ]  A$  + 

S“TR  + 

[  (pUv)i+1/2j  -  (pUv)i_1/2 j  ]  All  +  [  (pVv)iJ+1/2  -  (pVv)i j.1/2  ]  A£,  = 

t  (r9i  lfr  )i+i/2 j  -  (T<li  i?r  )i-i/2 j  1  Ar|  +  (4. 2) 

ds  ds 

[  (rq22^  )ij+l/2  -  (Tq22^  )ij-l/2 1  + 
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+  sv(tn) 

where  the  coordinate  transformation  employed  is  such  that  A£  =  At]  =  1. 

The  discrete  form  of  the  equations  governing  u  and  v  requires  knowledge  of  fluxes, 
exchange  coefficients,  and  dependent  variables  at  the  cell  faces.  Values  of  the 
transformation  metrics  are  also  required  at  the  cell  faces.  In  the  present  solution  procedure, 
some  of  the  transformation  metrics  are  considered  as  being  stored  at  the  cell  faces  and 
therefore  no  special  treatment  is  required.  For  the  determination  of  exchange  coefficients, 
fluxes,  and  dependent  variables,  various  interpolation  techniques  have  been  utilized. 

In  calculating  the  exchange  coefficient  at  a  particular  cell  face,  a  harmonic  mean  has 
been  used  (Patankar,  1980).  The  exchange  coefficient  at  the  i+1/2  face  (see  Figure  3)  is 
thus  determined  as  follows: 


■  i+l/2j 


r:5r 


1  i+lj 


fxi+l/2j  ^i+lj  +  (l'fxi+l/2j)  Fjj 


(4.  3) 


where  fxj+j/2j  is  the  £  coordinate  interpolation  factor  at  the  i+1/2  face  of  the  finite  control 
volume.  The  harmonic  mean  has  been  used  instead  of  linear  interpolation  because  of  the 
former's  ability  to  handle  abrupt  changes  in  the  distribution  of  the  exchange  coefficient. 
Such  a  characteristic  was  found  by  Nallasamy  and  Chen,  1985  to  be  crucial  in  obtaining 
physically  realistic  distributions  of  some  flow  field  variables. 

The  value  of  the  Cartesian  velocity  components  at  a  cell  face  is  obtained  by  assuming  a 
piecewise  linear  distribution  in  the  physical  space.  Such  a  procedure  is  consistent  with  a 
central  differencing  scheme.  For  large  values  of  cell  Peclet  number  (Pe  >  2)  a  central 
differencing  scheme  becomes  unstable  and  thus  the  cell  face  value  is  taken  to  be  equal  to  the 
velocity  of  the  upstream  neighbor  value.  The  use  of  central  differencing  for  low  Peclet 
number  and  upwinding  for  large  Peclet  number  is  known  as  hybrid  differencing  (Patankar, 
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1980).  Various  assumptions  can  be  made  as  regards  the  piecewise  distribution  of  the  flow 
variables.  Numerous  studies  have  been  made  to  ascertain  the  efficiency  of  some  of  these 
differencing  schemes  (Shyy  and  Correa,  1985;  Pollard  and  Siu,  1982).  Some  of  the  more 
popular  are  skew  upwind  differencing  (Raithby,  1976),  quadratic  upstream  interpolation 
for  convective  kinematics,  QUICK  (Leonard,  1979),  and  second  order  upwind 
differencing. 

The  hybrid  differencing  scheme  is  formally  only  first  order  accurate  and  is  thus 
questionable  as  regards  grid  independent  solutions.  More  significant  than  the  formal  order 
of  accuracy  is  the  false  diffusion  which  upwinding  introduces  into  the  calculation 
(Patankar,  1980).  However,  from  a  practical  point  of  view,  the  computational  efficiency  of 
a  solution  procedure  is  more  important  than  accuracy  alone  (Fletcher,  1988a).  The 
accuracy  of  even  a  formally  low  order  solution  may  be  improved  by  refining  the  grid. 
Also,  procedures  such  as  Richardson  extrapolation  can  be  employed  to  improve  solution 
accuracy  (Fletcher,  1988a).  The  hybrid  differencing  scheme  has  been  emplemented  in  the 
present  solution  procedure  because  it  is  unconditionally  stable. 

The  values  of  the  fluxes  at  the  cell  faces  are  not  determined  by  a  simple  linear 
interpolation.  In  the  context  of  a  collocated  grid  for  velocities  and  pressure,  such  a  practice 
would  invariably  lead  to  a  checkerboard  split  in  the  pressure  field  (Patankar,  1980).  The 
proper  treatment  of  the  cell  face  fluxes  is  crucial  for  the  development  of  a  convergent  and 
accurate  numerical  procedure.  The  cell  face  fluxes  are  ultimately  written  in  terms  of  the 
discretized  momentum  equations.  Therefore,  a  detailed  exposition  of  the  interpolation 
practice  used  for  cell  face  fluxes  is  deferred  until  the  final  forms  of  the  discretized 
momentum  equations  are  obtained. 

Given  the  interpolation  practices  for  the  various  cell  face  quantities,  the  discretized 
equations  for  u  and  v  are  given  as  follows: 
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Ap  uy  —  Ag  Ui+ij  +  Aw  Ui_i  j  +  An  Uij+i  +  As  Uij.i  +  Sjr  +  SU(^,T|)  (4.4) 

Ap  vy  =  Ag  vi+y  +  Aw  vj.i  j  +  An  vy+i  +  As  vjj.i  +  +  Sv(^,tj)  (4.  5) 

Ae  =  amaxl(  0.0,  -(pU)i+1/2j,  -(pU)i+1/2j  fxi+i/2j  +  (qn)i+i/2j  ri+i/2j  )  (4-6) 
Aw  =  amaxl(  0.0,  (pU)i-l/2j>  (pU)i-l/2j  (l"fxi-l/2j)  +  (Qll)i-l/2j  ^i-l/2j  )  (4-7) 

An  =  amaxl(  0.0,  -(pV)ij+1/2,  -(pV)jj+1/2  fyij+i/2  +  (<l22)io+l/2  rij+l/2 )  (4-  8) 

As  =  amaxl(  0.0,  (pV)ij.1/2,  (pV)ij.i/2  (1  -  fyij-1/2)  +  (<l22)ij-l/2  rij-l/2 )  (4-  9) 

Ap  =  Ae  +  Aw  +  An  +  Aj  (4. 10) 

S"fl;,ll>  =  -(an||+a21|^  (4.11) 

Svft.Tl)  —  ^  »12  ^  +  «J2  |^  )  (4.12) 

S™  =  |(r<)  +  J(rq21|)  (413) 

where  amaxl  is  a  fortran  supplied  function  which  is  used  in  implementing  the  hybrid 
differencing  procedure.  It  is  important  to  note  here  that  the  coefficients  for  u  and  v  are 
identical.  The  coefficients  for  u  also  serve  as  the  coefficients  for  v.  This  occurs  since  a 


42 


collocated  grid  is  used.  If  a  staggered  were  used,  then  the  coefficients  of  u  and  v  would  be 
different.  Indeed,  the  coefficients  of  the  discretized  equations  for  all  flow  variables  (k,  e, 
h0,  etc.)  are  the  same.  Therefore,  the  same  routine  can  be  used  for  the  calculation  of  the 
coefficients  of  any  flow  variable. 

In  looking  at  the  coefficients  of  the  discretized  equations  for  u  and  v,  it  is  evident  that 
all  of  the  coefficients  are  positive.  This  is  necessary  for  the  numerical  procedure  to  remain 
stable.  Also,  both  the  convective  and  diffusive  fluxes  at  a  given  cell  face  are  not  ascribed  to 
a  particular  cell.  The  cell  face  fluxes  are  not  calculated  as  though  they  belong  to  any  given 
cell.  This  quality  ensures  that  inconsistencies  do  not  arise  in  the  discretized  equations  for 
adjacent  control  volumes.  The  flux  leaving  one  control  volume  is  assured  of  being  the  flux 
which  enters  the  neighboring  control  volume.  Finally,  the  central  coefficient  (Ap)  is  seen  to 
be  equal  to  the  sum  of  the  neighboring  coefficients: 

Ap  =  Ae  +  Aw  +  An  +  As  +  {  (pU)i+1/2j  -  (pU)i_1/2j  +  (pV)ij+1/2  -  (pV)ij.i/2  }  (4. 15) 

==>Ap  =  ZAnb  (4.16) 

where  the  term  enclosed  in  brackets  is  the  discretized  continuity  equation  and  is  equal  to 
zero.  This  condition  is  considered  necessary  for  situations  in  which  the  differential 
equation  is  satisfied  even  if  a  constant  were  to  be  added  to  the  dependent  variable 
(Patankar,  1980). 

When  solving  for  u  or  v  (or  any  flow  variable)  in  an  iterative  manner,  underrelaxation 
is  employed  in  order  to  maintain  stability  of  the  numerical  procedure.  If  underrelaxation  is 
not  employed,  the  scheme  can  become  unstable  since  the  coupled  nonlinear  equations  are 
treated  in  a  decoupled,  linear  fashion.  In  the  present  solution  procedure,  the 
underrelaxation  is  implemented  in  an  implicit  manner  (Patankar,  1980).  The  final  form  of 
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the  discretized  equations  for  u  and  v  as  used  in  the  present  solution  procedure  are  thus 
given  as  follows: 


Ap  ug  =  Ag  ui+1J  +  Aw  uM  j  +  An  Ui  J+1  +  As  uy_!  +  +  S«(^,T|)  (4.  17) 

Ap  vij  =  Ae  vi+l  j  +  Aw  vi-l  j  +  An  Vi  j+1  +  As  Vy_!  +  +  Sv(£,Tl)  (4.  18) 


Ae  =  amaxl(  0.0,  -(pU)i+1/2J,  -(pU)i+1/2J  fxi+1/2J  +  (qn)i+1/2J  ri+1/2j )  (4.  19) 
Aw  =  amaxl(  0.0,  (pU)M/2j,  (pU)i.1/2J  (l-fxi.1/2J)  +  (qn)^^  Tj.^j )  (4.20) 
An  =  amaxl(  0.0,  -(pV)y+1/2,  -(pV)y+1/2  fyiJ+1/2  +  (q22)y+i/2  riJ+1/2 )  (4.  21) 
As  =  amaxl(  0.0,  (pV)y_1/2,  (pV)y.1/2  (1  -  fyy_1/2)  +  (q22)ij.1^  ry.1/2 )  (4.  22) 


(4.23) 

(4.24) 

(4.  25) 

(4.  26) 

(4.27) 
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In  the  present  solution  procedure,  the  method  of  Rodi,  Majumdar,  and  Schonung,  1987 
has  been  used  to  determine  the  cell  face  values  of  the  convective  flux.  Various  methods 
exist  for  the  treatment  of  collocated  or  non-staggered  grids  (Rhie  and  Chow,  1983;  Thiart, 
1990;  Acharya  and  Moukalled,  1989;  Miller  and  Schmidt,  1988).  The  present  method  is 
similar  to  that  employed  by  Rhie  and  Chow.  It  differs  from  their  method  in  that  it  is  not 
relaxation  dependent  (Majumdar,  1988).  Also,  Acharya  and  Moukalled,  1989  found  that 
the  method  proposed  by  Rhie  and  Chow  did  not  satisfy  integral  mass  conservation  because 
of  the  flux  correction  which  is  employed  in  that  method.  The  primary  concept  of  the 
method  proposed  by  Rhie  and  Chow  was  that  pressure  oscillations  should  not  be  allowed 
to  occur  as  evidenced  by  physical  reasoning.  This  feature  is  retained  in  the  present  solution 
procedure. 

Following  Majumdar  (Majumdar,  1988)  the  values  of  the  Cartesian  velocities  at  the  cell 
faces  are  written  in  terms  of  the  interpolated  discretized  equations  for  the  adjacent  cell  center 
velocities.  For  example,  if  a  staggered  grid  were  being  used,  the  values  of  ue  and  ve  would 
be  written  as  follows: 


ue  = 


_  (^^nbunh)e  \  \ 


\  ^  J  Kl 


\  ’a"  af  +  821  dn  >c  |a"~  S"®’T1)  [  +  < 1  •  “  >  “e  d  IA  28> 


ve  = 


'(^AnbVnb)e 


1  ^  J 


J. 


■ ,  ap ,  ap ,  , 

^(ai2aTa2We^ 


V 


S.  +  (l-©)v°ld  (4.29) 


where  SU(^,T])  and  Sv(^,t|)  obviously  do  not  contain  the  source  term  due  to  pressure  or 
relaxation.  Now,  a  staggered  is  not  being  used  and  so  quantities  such  as  Ap  and  ZAnbun^ 
are  not  known  at  the  cell  faces.  To  determine  these  quantities,  linear  interpolation  in  the 
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physical  space  is  used.  In  particular,  only  those  terms  surrounded  by  braces  are 
interpolated.  The  pressure  terms  and  relaxation  terms  are  not  interpolated.  The  cell  face 
contravariant  velocity  in  terms  of  ue  and  ve  is  as  follows: 


Ue=  uean  +  vea12 


(4.  30) 


where  aj  j  and  a12  have  been  calculated  at  the  cell  face  (east  in  this  case).  With  the 
interpolation  of  the  various  terms  in  ue  and  ve  included,  the  cell  face  contravariant  velocity 
is  thus  expressed  in  terms  of  cell  center  quantities: 


ue  =  all 


+  al2 


■^^nbunb  +  SU(^,T|)  ^ 

-an 

1  1 

K 

p  j 

AP 

J 

zAnbvnh  +  sv(tn)  ] 

-  a12 

f~T  \ 

*  J 

'i 

k  y 

(all^  +  a21  f^)e  +  an(l-co)ue0,d 

an 


,  ap  _ ap , 

(a127T  +  a22  “)e 
OS  oq 


+  a12  (1  -  co)  veold 


(4.31) 


where  the  overbar  indicates  that  a  linear  interpolation  in  the  physical  space  has  been  used  to 
obtain  that  quantity.  If  the  terms  involving  ue  and  ve  are  grouped  together,  the  resulting 
expression  for  Ue  shows  that  it  will  be  independent  of  the  relaxation  factor  when  a 
converged  solution  is  obtained: 


Ue  =  co  Uenew  +  (1  -  co )  Ueold  (4.32) 

The  expressions  for  Uw,  v„,  and  Vs  follow  a  similar  logic.  This  procedure  for  determining 
cell  face  velocities  has  been  referred  to  as  momentum  interpolation  (Majumdar,  1988). 
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The  cell  face  value  of  the  contravariant  velocity  has  been  represented  in  terms  of  cell 
center  values  but  nothing  has  been  mentioned  as  regards  density.  The  value  of  the  cell  face 
density  is  obtained  by  a  strict  upwinding  procedure  (Issa  and  Lockwood,  1977).  In  the 
present  method  then,  the  cell  face  flux  is  not  treated  as  one  entity.  The  cell  face  flux  is 
determined  by  the  product  of  the  upwind  density  and  the  contravariant  velocity  as  obtained 
from  momentum  interpolation. 

4.2.  Continuity  Equation 

In  the  present  procedure,  mass  conservation  is  satisfied  by  applying  corrections  to  the 
velocities  (contravariant  and  Cartesian).  The  discretized  form  of  the  mass  conservation 
statement  serves  as  the  basis  for  the  development  of  a  suitable  correction  equation. 


<pU>i+i/2 j  -  (pUViflj  +  (pV)ij+1/2  -  (pV)ij.  1«  =  0  (4.  33) 

The  velocities,  density,  and  pressure  are  written  as  the  sum  of  their  estimated  values  and  a 
suitable  correction. 


p  =  p*  +  p' 

(4.34) 

p  =  p*  +  p 

(4. 35) 

U  =  U*+U’ 

(4.  36) 

V  =  v*  +  V' 

(4.  37) 
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Upon  substitution  into  the  discretized  continuity  equation,  the  following  is  obtained: 


(P*U‘  +  P’U*)i+1/2J  -  (p*U(  +  p'U*)^  +  (p*V  +  p’V*)ij+i/2  -  (P*V'  +  p'V*)i,j.i/2 
=  -  [  (P*U*)i+1/2j  -  (p*U*)i_1/2J  +  (p*V*)iJ+1/2  -  (p*V*)i(j.1/2  ]  (4.  38) 

It  has  been  assumed  that  the  nonlinear  compressibility  terms  (p'U',  p'V')  are  negligible. 
In  the  work  of  Shyy  and  Braaten,  1988  these  terms  were  included  in  the  pressure 
correction  equation  as  an  explicit  source.  It  was  found  that  inclusion  of  these  terms  was 
useful  in  preventing  early  divergence  of  the  numerical  procedure  for  the  calculation  of 
supersonic  flows.  The  present  method  is  envisioned  for  use  in  calculating  subsonic 
compressible  flows  and  so  the  nonlinear  terms  were  not  considered  necessary.  When 
incompressible  or  low  mach  number  flows  are  calculated,  the  density  corrections  are 
removed  altogether. 

The  right  hand  side  of  the  above  correction  equation  represents  a  mass  source.  When  a 
converged  solution  is  obtained,  the  velocities  as  calculated  from  the  momentum  equations, 
together  with  the  density,  will  be  such  that  this  mass  source  is  zero.  Since  corrections  are 
not  applied  to  the  converged  result  from  the  momentum  equations,  the  correction  equation 
does  not  effect  the  formal  accuracy  of  the  procedure.  The  correction  equation  does 
however  affect  the  convergence  rate  of  the  overall  numerical  procedure.  It  is  therefore 
vitally  important  to  utilize  a  correction  equation  which  is  convergent.  The  issue  of  a 
convergent  correction  equation  arises  when  considering  the  form  taken  by  the  corrections  in 
a  non-orthogonal  coordinate  system. 

Each  of  tiie  velocity  and  density  corrections  is  written  in  terms  of  a  pressure  correction. 
In  this  way  the  discretized  continuity  equation  becomes  a  pressure  correction  equation.  The 
mass  source  in  the  equation  serves  to  drive  these  pressure  corrections  and  thus  satisfy 
continuity. 
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The  density  correction  as  a  function  of  pressure  is  obtained  from  the  isentropic  speed  of 
sound  for  an  ideal  gas: 


dP 

aP 


=  yRT 


therefore, 


P’  = 


(4.  39) 


(4.  40) 


Other  researchers  (Shyy  and  Braaten,  1988;  Issa  and  Lockwood,  1977;  Demirdzic,  et.al., 
1990)  have  obtained  the  expression  for  p'  from  the  ideal  gas  law,  hence 


P'  = 


(4.41) 


The  use  of  the  isentropic  sound  speed  in  expressing  p'  is  thus  seen  to  be  merely  an 
underrelaxed  version  of  the  expression  obtained  from  the  ideal  gas  law.  The  specific  heat 
ratio  (y)  ensures  that  p'  as  obtained  from  the  sound  speed  definition  will  be  less  than  the  p’ 
obtained  from  the  ideal  gas  law. 

The  corrections  to  the  contravariant  velocities  (U'  and  V')  are  obtained  from  the 
definitions  of  the  contravariant  velocities  and  from  the  momentum  equations. 


U'  =  au  u’  +  &12  v' 

(4.42) 

V'  =  a2i  u’  +  a22  V 

(4.  43) 
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As  an  example,  consider  the  form  of  the  correction  for  the  contravariant  velocity  at  the  east 
face  of  a  control  volume  (U'i+i/2j)- 

U’i+l/2j  =  al  1  u'i+l/2j  +  a12  v'i+l/2j  (4.  44) 

The  transformation  metrics  ajj  and  aj2  are  stored  at  the  cell  faces  corresponding  to  lines  of 
constant  In  this  way,  no  averaging  need  be  performed  in  order  to  express  these 
quantities  at  the  east  or  west  face  of  a  finite  control  volume.  The  values  of  u'  and  v’  are 
determined  from  the  momentum  equations  in  a  manner  consistent  with  the  SIMPLEC 
method  suggested  by  VanDoormaal  and  Raithby,  1984.  Consider  the  form  of  the 
discretized  equations  for  the  estimated  velocities  (u*  and  v*)  at  the  east  face  of  a  control 
volume: 


[a;].  u*e  =  [£AnbU.nb]e -<anfpa21|-\ 

+  [Sra]e+ Met 1 (4.45) 

[ApV]eV*e  =  [lAnbv.nb]e.(a12^  +  a22|i)e 

+  [STRle+MeO-^Jvf1  (4.46) 

Since  a  collocated  grid  is  being  used,  the  coefficients  (Ap,  Anb)  are  not  defined  at  cell  faces. 
These  coefficients  are  obtained  by  linear  interpolation  of  the  cell  center  quantities  as 
mentioned  previously.  If  these  equations  for  u*e  and  v*e  are  subtracted  from  the 
corresponding  equations  for  ue  and  ve  (eqns.  4.28  and  4.29),  then  the  following  equations 
for  u’e  and  v’e  are  obtained. 
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(4. 47) 


[Aj]e  u'e  =  [XAnbu*nb]e  -  ( an||  +  a21|^ )e 

[Ap]e  v'e  =  [£Anbv'nb]e  -  ( a12|^  +  a22|^  )e  (4.  48) 

In  accordance  with  the  SIMPLEC  method,  the  quantity  (EA“b)u'e  is  subtracted  from  both 
sides  of  the  equation  for  u'e  and  (LA^b)v'e  is  subtracted  from  both  sides  of  the  equation  for 

v'e.  The  final  expressions  for  u'e  and  v’e  are  obtained  by  disregarding  the  terms  involving 
the  neighbor  coefficients  (VanDoormaal  and  Raithby,  1984).  The  form  of  the  equation  for 
u'e  and  v'e  is  therefore, 


ue  =  - 


(Ap  -  LAnb) 


3P’  8P'  \ 

aT'321^) 


ve  =  ‘ 


(A, 


= — -  f 

-  2Anb)  l 


ap'  ap*  n 
_ai2  aTa22^) 


(4. 49) 


(4.50) 


where  the  coefficients  (Ap,  Anb)  are  obtained  from  a  linear  interpolation  of  cell  center 
values.  It  is  to  be  realized  that  the  removal  of  the  neighbor  coefficient  terms  from  the 
equations  for  u'e  and  v'e  is  a  legitimate  procedure.  In  the  converged  state,  all  pressure 
corrections  will  be  zero  since  the  mass  source  is  zero.  The  converged  solution  is  then 
uninfluenced  by  the  form  of  the  pressure  correction  equation  since  it  is  not  used  in  the  final 
iteration  (Patankar,  1980).  The  form  of  the  correction  which  is  applied  to  U'e  is  obtained 
then  by  substitution  of  u'e  and  v'e  into  equation  4.44: 

u'e = L  L  r\ ( ai‘2 + ai22  >e  ( p'y  ■  p'i+iJ } 

^  (Ap  -  £Anb)  j 
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(4. 51) 


(At 


^Anb)  j 


(all  a21  +a12  a22  )e 


(-5> 


where  the  overbar  denotes  an  interpolated  quantity.  The  expressions  for  U'w,  V'n,  and  V's 
are  obtained  in  an  analogous  manner. 

Note  in  the  above  expression  for  U'e  that  two  pressure  gradients  are  present,  one  in  the 
coordinate  direction  and  one  in  the  T)  direction.  The  pressure  gradient  in  the  rj  direction  in 
the  above  expression  for  U'e  is  zero  for  an  orthogonal  grid.  It  is  only  present  when  non- 
orthogonal  grids  are  used.  This  fact  has  led  some  researchers  (Rhie  and  Chow,  1983; 
Acharya  and  Moukalled,  1989)  to  neglect  these  terms,  thus  assuming  near  orthogonality. 
As  noted  above,  the  assumptions  made  in  deriving  the  pressure  correction  equation  do  not 
effect  the  final  solution,  only  the  rate  of  convergence  of  the  solution  procedure.  Neglecting 
these  non-orthogonal  terms  is  thus  a  legitimate  procedure.  However,  Peric,  1990  has 
found  that  the  applicability  of  a  solution  procedure  which  incorporates  the  near  orthogonal 
simplification  is  severely  restricted  when  the  grid  becomes  non-orthogonal.  For  this 
reason,  the  present  solution  procedure  retains  these  non-orthogonal  terms  in  the  pressure 
correction  equation.  The  inclusion  of  the  non-orthogonal  terms  results  in  a  nine  point 
computational  molecule  for  the  pressure  correction  equation.  Rather  than  use  a  more 
complex  solver,  the  non-orthogonal  terms  have  been  included  as  explicit  source  terms,  thus 
allowing  the  use  of  only  a  five  point  computational  molecule.  Therefore,  the  final  form  of 
the  pressure  correction  equation  as  used  in  the  present  procedure  is  as  follows: 

Aj  P'i j  =  Ae  j  +  Aw  pi-l j  +  AS  J-yi  +  A?  P'ij-1  +  SM  +  SNO  (4.  52) 


( Ap-£A„b  ) 


i+l/2j 


2  2 

all+  a12  )i+l/2j  pi+l/2j 
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-  UVlflj  fxi+i/2j  ^)+l  j  (4. 53) 

A“  =  [<Ap-Lnb)]  M«  '  S'1+  a‘2  )M/2J  PM« 

+  u*i-l/2j  ( 1  -  fxi-l/2j )  ;Yi  j  (4-  54) 

A"  =  [(  Ap-IAnb  )]  *0+1/2  (  $2  \ j+1/2  Pi  j+1/2 

-  V*ij+i/2  fyij+1/2  (^j+i  (4. 55) 

A"  =  [(  Ap-Lnb  )]  >0-1/2  (  ^1+  ^2  >iJ-l/2  PiJ-I/2 

+  v*ij-l/2  ( 1  -  fyij-i/2 )  (4.  56) 

AP  =  [(  Ap-Linb  )]  l+'flo  ( a"+  *“  >W«  Pi+1«J 
+  =(  Ap-ZAnb  )]  >-W0  ( “>1+  a'2  )>->«*i  P>-1«0 
_(  Ap-EAnb  )  ij+1/2  ( ^l+  a 22  )iJ+1/2  Pi>j+1/2 


i  I  n 

_(  Ap-LAnb  )_  io-1/2 ( a21+  322  )lJ'1/2  PiJ’1/2 
+  U*i+1/2j  ( 1  -  fxi+1/2j  )  -  U*i.^j  fxi.1/2j 

+  V**j+l/2  <  1  *  *1*1/2  )  (-^j  i j  *  V*i j-i/2  fyij-l/2  j  (4.  57) 


SM  =  -  [  (P*U*)i+l/2j  -  (p*U*)j.1/2j  +  (p*V*)iJ+i/2  -  (p*V*)jj.1/2  ]  (4.  58) 


SN0  = 
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+ 


L  J 


(  an  a2l  +  ai2  a22  )i-l/2j  (^")-l/2j 
(  a21  all  +  a22  a12  )y+l/2  (^)j+l/2 
(  a21  all  +  a22  a12  >i j-l/2 


(4.  59) 


4.3.  General  Scalar  Equation 

After  the  momentum  and  continuity  equations  are  solved,  the  solution  to  the  scalar 
variables  (k,e,YFU,  etc.)  must  be  obtained.  Since  the  scalar  variables  are  stored  at  the  cell 
centers  along  with  the  Cartesian  velocity  components,  the  discretized  equation  for  u  and  v 
also  serves  as  the  discrete  equation  governing  the  distribution  of  any  flow  variable.  Only 
the  diffusive  exchange  coefficient  and  the  source  term  will  change  for  the  solution  to  any 
particular  flow  variable.  The  exchange  coefficients  and  source  terms  for  the  various  scalar 
quantities  used  in  the  present  numerical  algorithm  are  given  in  Table  1. 

In  solving  for  the  scalar  variables,  it  was  necessary  to  employ  source  term  linearization 
(Patankar,  1980).  All  portions  of  the  source  terms  which  were  positive  were  retained  on 
the  right  hand  side  of  the  discrete  equation.  All  negative  quantities  in  the  source  term  were 
taken  to  the  left  hand  side  of  the  discrete  equation  and  incorporated  into  the  central 
coefficient  (Ap).  This  practice  was  found  to  be  vitally  important  in  ensuring  that  turbulent 
kinetic  energy  and  turbulent  dissipation  rate  not  become  negative. 

4.4.  Solution  Procedure  Summary 

As  mentioned  previously,  the  solution  procedure  employed  involves  a  sequential 
solution  of  various  flow  variables.  Figure  4  presents  the  steps  taken  in  obtaining  a 
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converged  solution  to  a  given  flow  problem.  The  procedure  involves  solving  for  the 
Cartesian  velocity  components  u  and  v.  The  pressure  correction  equation  is  then  solved 
with  a  subsequent  correction  of  the  relevant  flow  variables  (u,  v,  U,  V,  and  P).  The  scalar 
variables  are  solved  next.  The  order  in  which  the  scalars  are  solved  is  as  follows: 

*  Swirl  Velocity  Component  (w) 

*  Turbulent  Kinetic  Energy  (k) 

*  Turbulent  Dissipation  Rate  (e) 

*  Mixture  Fraction  (f) 

*  Fuel  Mass  Fraction  (Yfu) 

*  Concentration  Fluctuation  (g) 

*  Stagnation  Enthalpy  (h0) 

The  final  step  in  the  iterative  sequence  involves  an  update  of  thermodynamic  and  transport 
properties. 

The  solution  of  any  particular  flow  variable  involves  the  solution  of  a  set  linear 
algebraic  equations.  An  alternating  direction  implicit  (ADI)  method  has  been  used  for  this 
purpose.  The  ADI  uses  a  tridiagonal  matrix  algorithm  (TDMA)  to  obtain  a  direct  solution 
to  all  variables  at  a  given  £,  coordinate  location  (Anderson,  et.  al.,  1984).  After  sweeps  are 
performed  in  the  2;  direction,  a  similar  procedure  is  performed  in  the  Tj  direction.  A  certain 
number  of  sweeps  through  the  ADI  was  performed  for  each  flow  variable.  In  general,  5 
sweeps  were  used  for  the  velocity  components  and  scalars  while  10-20  sweeps  were  used 
for  the  pressure  correction. 

In  testing  the  present  solution  procedure,  it  was  found  beneficial  to  apply  integral  mass 
adjustments  to  the  cell  face  fluxes  for  some  flow  problems.  Two  forms  of  implementation 
of  the  integral  mass  adjustment  were  tested.  One  procedure  follows  from  the  work  of 
Patankar  and  Spalding,  1972  and  involves  adding  a  constant  correction  to  all  the 
contravariant  velocities  in  a  given  plane  so  as  to  satisfy  continuity.  The  pressure  was  also 
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corrected  so  as  to  reflect  the  change  in  the  mass  flux.  The  second  procedure  involved 
scaling  the  fluxes  at  a  given  downstream  plane  so  that  the  total  mass  flow  rate  through  the 
upstream  and  downstream  planes  was  identical. 


’ny+1 

"ny+1 

S  (PU)MJ 

+  (PV)i,l  +  (pV)j>ny+l  —  0C 

X(pU)ij 

L  j=2  J 

L  j=2  J 

(4.  60) 


The  factor  a  was  calculated  at  a  given  plane  based  on  the  estimated  mass  flux  through  that 
plane.  All  fluxes  were  then  multiplied  or  scaled  by  the  factor  a  so  as  to  satisfy  integral 
mass  conservation.  A  corresponding  pressure  correction  was  subsequently  added  to  all  the 
pressures  on  the  downstream  side  of  the  given  plane. 


(4.  61) 


By  applying  the  integral  mass  adjustments  to  all  the  planes  perpendicular  to  the  main  flow 
direction,  overall  mass  conservation  is  satisfied. 

Through  tests  of  the  computational  procedure,  it  was  found  that  the  second  method  of 
integral  mass  adjustment  worked  better.  This  is  presumably  due  to  the  fact  that  the  addition 
of  a  correction  to  a  given  flux  may  cause  it  to  change  sign.  Such  an  occurence  would  not 
be  valid  in  a  region  of  flow  reversal.  Further  observations  as  regards  integral  mass 
adjustments  and  on  the  solution  procedure  as  a  whole  are  presented  in  the  following 
chapter. 
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CHAPTER  5:  EVALUATION  OF  NUMERICAL  METHOD 


It  is  essential  to  know  the  applicability  and  efficiency  of  any  flowfield  simulation.  The 
applicability  is  assessed  by  the  accuracy  of  predictions  made  with  the  simulation.  The 
efficiency  is  gauged  in  part  by  the  convergence  rate  and  how  it  is  affected  by  such  factors 
as  cell  aspect  ratio,  grid  density,  grid  distribution,  and  flow  Reynolds  number.  The 
evaluation  of  the  present  simulation  of  turbulent  reacting  flow  is  made  by  comparing  its 
predictions  with  experimental  results  found  in  the  literature.  Six  flow  configurations  have 
been  considered  in  evaluating  the  numerical  simulation:  developing  turbulent  flow  in  a 
straight  pipe,  turbulent  flow  in  an  axisymmetric  sudden  expansion,  laminar  diffusion 
flame,  turbulent  diffusion  flame,  turbulent  premixed  flame,  and  turbulent  swirl  flow  in  an 
axisymmetric  sudden  expansion.  The  flowfield  predictions  are  compared  with 
experimental  results.  Conclusions  are  made  concerning  the  behavior  of  the  numerical 
simulation  so  as  to  assess  the  efficiency  and  applicability  of  simulating  turbulent  reacting 
flows  with  the  solution  procedure. 

5.1.  Turbulent  Pipe  Flow 

The  developing  turbulent  flow  in  a  pipe  of  circular  cross  section  has  been  calculated 
using  the  numerical  algorithm.  The  results  of  the  computations  have  been  compared  with 
the  experiments  of  Richman  and  Azad  (Richman  and  Azad,  1973)  and  C.  J.  Lawn 
(Martinuzzi  and  Pollard,  1989).  For  the  purpose  of  calculating  the  flowfield,  the  pipe  was 
taken  to  be  40  units  long  with  a  radius  of  0.5  units. 

The  calculation  was  made  by  assuming  the  existence  of  a  uniform  velocity  profile  at  the 
pipe  inlet.  The  inlet  conditions  for  turbulent  kinetic  energy  and  dissipation  rate  were  also 
given  as  uniform  at  values  of  0.005  U2  and  ( Cpk3/2  /  0.03R )  respectively  (Martinuzzi  and 
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Pollard,  1989).  Since  the  flow  was  axisymmetric,  the  computational  domain  consisted  of 
an  axis  of  symmetry  at  the  lower  boundary  and  a  wall  at  the  upper  boundary.  A  fully 
developed  flow  for  all  flow  variables  was  prescribed  at  the  outflow  plane.  The  use  of  the 
hybrid  differencing  scheme  however,  does  not  allow  the  information  at  the  outflow  plane  to 
propagate  into  the  inner  flow  domain.  The  outflow  boundary  condition  is  thus  one  of 
parabolic  flow  (Patankar,  1980). 

The  comparison  of  experimental  and  numerically  obtained  mean  velocity  profiles  for  a 
Reynolds  number  of  300,000  is  given  in  Figures  5  through  9.  The  mean  streamwise 
velocity  component  (u)  has  been  non-dimensionalized  with  the  bulk  velocity  (U)  and  has 
been  plotted  as  a  function  of  the  radius  ratio  (R  is  the  pipe  radius).  The  calculation  has 
been  performed  using  an  80  x  40  grid.  The  cell  size  in  the  axial  direction  up  to  x/D  of  2 
was  given  as  0.35  while  the  cell  size  thereafter  was  equal  to  0.5.  Therefore,  only  slight 
concentration  of  the  mesh  in  the  inlet  region  was  utilized.  It  is  observed  that  the  calculation 
overpredicts  experiment  in  the  region  close  to  the  wall.  The  discrepancy  is  seemingly  due 
to  insufficient  grid  resolution  in  the  axial  direction.  The  calculated  velocity  profiles  do 
progress  toward  the  fully  developed  condition.  The  distance  over  which  this  development 
occurs  (entrance  length)  is  overpredicted  with  the  use  of  an  80  x  40  grid.  The  comparison 
reveals  that  the  discrepancy  diminishes  as  the  flow  develops  (Figure  9).  The  use  of  80 
cells  in  the  axial  direction  is  not  entirely  sufficient  to  capture  the  steeper  gradients  in  mean 
flow  velocity  which  occur  in  the  inlet  due  to  the  merging  boundary  layers.  A  non-uniform 
grid  with  increased  mesh  concentration  in  the  inlet  region  may  prove  more  appropriate. 

The  streamwise  velocity  development  at  a  given  radius  is  seen  in  Figures  10  through 
12.  These  results  are  consistent  with  the  high  Reynolds  number  k  -  e  model  prediction 
given  by  Martinuzzi  and  Pollard,  1989.  The  prediction  is  observed  to  be  better  in  the  core 
than  in  the  near  wall  region. 
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The  dimensionless  turbulent  kinetic  energy  comparison  given  in  Figure  13  is  quite  good 
(u*  is  the  friction  velocity).  The  experimental  data  are  for  fully  developed  turbulent  pipe 
flow  while  the  numerical  results  are  for  x/D  equal  to  40.  The  mild  difference  observed  may 
be  due  to  the  fact  that  x/D  equal  to  40  may  not  be  fully  developed  for  a  Reynolds  number  of 
380,000  and  hence  the  numerical  results  presented  are  not  for  fully  developed  flow.  These 
results  are  comparable  to  those  obtained  by  Martinuzzi  and  Pollard,  1989  for  a  high 
Reynolds  number  k-e  model. 

These  results  reveal  that  the  numerical  algorithm  presented  can  be  used  to  predict  not 
just  the  qualitative  features  of  turbulent  pipe  flow  but  the  quantitative  features  as  well.  The 
quantitative  accuracy  however,  is  dependent  upon  a  sufficient  resolution  of  the 
computational  grid. 

5.2.  Turbulent  Flow  in  an  Axisymmetric  Sudden  Expansion 

The  experiments  of  Craig,  et.  al.,  1984  are  used  for  comparative  purposes.  The 
geometry  considered  is  that  of  an  axisymmetric  sudden  expansion  of  area  ratio  3.57.  The 
flow  Reynolds  number  at  the  inlet  was  86,293.  Experimental  data  was  obtained  using  a 
laser  doppler  velocimeter  with  consideration  given  to  the  velocity  biasing  which  occurs 
(Craig,  et.  al.,  1984).  Further  details  of  the  experiment  and  the  procedures  followed  are 
given  by  Craig,  et.  al.,  1984. 

The  simulation  of  the  flow  in  the  sudden  expansion  was  considered  to  be  axisymmetric. 
This  assumption  is  supported  by  the  experimental  findings  (Craig,  et.  al.,  1984).  It  was 
further  assumed  that  the  flow  was  time  mean  steady.  The  lower  boundary  of  the 
computational  grid  (Figure  51)  is  thus  specified  as  a  symmetry  boundary.  The  inlet 
velocity  was  specified  as  being  uniformly  distributed  at  27  m/s.  The  inlet  profiles  of 
turbulent  kinetic  energy  and  turbulent  dissipation  rate  were  likewise  prescribed  as  uniform. 
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The  turbulent  kinetic  energy  was  given  as  0.0009  and  the  dissipation  rate  at  1.2(10-5) 

in  accordance  with  the  simulation  of  Vanka,  1988.  The  outflow  plane  was  placed  at  a 
distance  of  16H  (H  is  the  step  height)  from  the  step.  The  flow  at  the  outlet  was  considered 
to  be  fully  developed  as  described  in  a  previous  chapter. 

The  comparison  of  experimental  and  calculated  profiles  for  the  mean  velocity  and 
turbulent  kinetic  energy  reveal  a  good  qualitative  agreement  (Figures  14  to  21).  The  most 
notable  disparities  in  experimental  and  calculated  profiles  occur  in  two  regions,  in  the 
recirculation  region  behind  the  expansion  and  near  the  exit  of  the  combustor. 

The  mean  velocity  profile  at  x/H  equal  to  13.5  is  overpredicted  for  the  use  of  40  x  20 
and  50  x  30  grids  (Figure  17).  The  calculation  on  a  60  x  40  grid  resulted  in  a  similar  result 
even  with  further  grid  refinement  in  the  vicinity  of  the  combustor  exit. 

Quantitative  disagreement  is  also  observed  in  the  flow  recirculation  region  behind  the 
step.  In  fact,  the  qualitative  agreement  is  erroneous.  The  experiment  shows  an  increase  in 
reverse  flow  velocity  behind  the  step  (Figure  14).  The  reverse  flow  velocity  subsequently 
decreases  as  the  edge  of  the  shear  layer  is  reached.  The  calculation  on  all  three  grids 
reveals  a  monotonic  decrease  in  the  magnitude  of  the  reverse  flow  velocity.  This 
observation  was  also  made  by  Obi,  et.  al.,  1991  when  a  k-£  turbulence  model  was  used. 
They  did  not  observe  this  disagreement  for  use  of  a  second  moment  closure. 

The  profiles  of  turbulent  kinetic  energy  all  reveal  that  the  simulation  underpredicts  this 
quantity  (Figure  18  to  21).  The  predicted  profiles  of  turbulent  kinetic  energy  reveal  that 
radial  diffusion  of  k  does  not  occur  soon  enough  along  the  axis  of  the  duct.  Even  at  x/H 
equal  to  7.8  (Figure  20),  the  centerline  turbulent  kinetic  energy  is  still  approximately  zero. 
The  experimental  results  are  seen  to  exhibit  more  diffusive  effects  than  the  calculation. 

Figure  22  gives  the  convergence  histories  corresponding  to  the  three  grids  used  in  the 
simulation.  The  solutions  were  assumed  to  be  converged  when  the  mass  residual  was 
reduced  by  5  orders  of  magnitude.  As  expected,  finer  meshes  result  in  more  iterations  to 
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reach  convergence.  These  calculations  were  made  without  considering  the  non-orthogonal 
terms  in  the  pressure  correction  equation.  Inclusion  of  these  terms  may  result  in  a  slightly 
better  convergence.  Also,  the  harmonic  mean  determination  for  cell  face  values  of 
exchange  coefficient  was  not  used.  In  the  context  of  a  non-orthogonal  grid,  it  led  to 
convergence  difficulties  and  in  some  instances  divergence. 

The  convergence  rates  given  are  for  varying  relaxation  factors.  The  solution  on  finer 
grids  required  more  relaxation  in  general.  The  relaxation  factors  for  velocity,  pressure,  and 
turbulence  quantities  were  0.7, 1.0,  and  0.8  for  use  of  a  40  x  20  grid,  0.7, 0.9,  and  0.7  for 
the  50  x  30  grid,  and  0.6, 0.5, 0.6  for  the  calculation  on  the  60  x  40  grid. 

5.3.  Laminar  Diffusion  Flame 

The  usefulness  of  the  one  step  infinite  rate  chemistry  assumption  for  modelling  laminar 
diffusion  flames  was  assessed  by  simulating  the  experiment  of  Mitchell,  et.  al.,  1980.  The 
flow  configuration  involved  a  confined  laminar  methane  diffusion  flame.  The  overall 
equivalence  ratio  was  approximately  0.3.  Radial  profiles  of  temperature,  velocity,  and 
specie  mole  fraction  were  obtained  at  various  axial  locations.  Comparison  with  these 
results  gives  an  indication  of  the  accuracy  of  the  present  numerical  simulation. 

The  geometry  for  the  given  problem  is  a  simple  pipe.  The  flow  was  assumed  to  be 
axisymmetric.  Inlet  velocity  profiles  for  both  the  fuel  and  the  air  were  taken  to  be  uniform. 
The  fuel  velocity  was  specified  at  4.5  cm/s  and  the  air  at  9.88  cm/s.  The  outflow  plane  was 
specified  at  13  tube  diameters  (13D)  from  the  inlet.  The  flow  was  assumed  to  be  fully 
developed  at  this  location. 

The  simulation  used  an  80  x  46  grid.  The  grid  was  refined  near  the  inlet  so  as  to  better 
resolve  the  physics  of  the  flame.  Also,  the  simulation  utilized  the  harmonic  mean  for 


61 


interpolating  the  laminar  viscosity.  A  laminar  Schmidt  number  of  unity  was  used  in  the 
simulation. 

Figures  23  through  25  reveal  a  qualitatively  good  agreement  with  experiment.  The 
calculated  temperatures  are  greater  that  experiment  due  to  radiation  and  finite  rate  effects. 
The  fuel  mole  fraction  profiles  (Figure  29  and  30)  show  that  the  flame  width  is  captured 
well.  Consideration  of  fuel  pyrolysis  in  the  numerical  simulation  would  improve  the 
prediction  of  fuel  mole  fraction.  The  fast  chemistry  model  tacidy  assumes  that  all  unreacted 
carbon  is  bound  up  in  the  fuel.  In  reality,  this  not  the  case  (Mitchell,  et.  al.,  1980). 
Figures  31a  and  b  give  color  contours  of  temperature  and  velocity.  The  flame  shape  as 
well  as  the  predicted  outer  flow  recirculation  are  readily  apparent 

The  predicted  velocity  profiles  all  agree  well  with  experiment  (Figures  26  to  28).  The 
outer  flow  recirculation  was  also  captured  by  the  simulation  (Figure  31b).  The  presence  of 
the  recirculation  region  was  the  reason  for  use  of  such  a  long  calculation  domain. 
Placement  of  the  outflow  boundary  in  the  region  of  the  recirculation  region  would  be  an  ill- 
prescribed  boundary  condition  and  would  most  likely  lead  to  divergence  of  the  numerical 
procedure. 

5.4.  Turbulent  Diffusion  Flame 

The  turbulent  diffusion  flame  studied  by  Lewis  and  Smoot,  1981  was  analyzed  using 
the  present  solution  procedure.  The  experiment  which  was  simulated  involved  the  issuing 
of  coaxial  jets  of  fuel  and  air  into  an  axisymmetric  sudden  expansion.  The  fuel  used  in  the 
experiment  was  a  mixture  of  methane  and  ethane.  The  momentum  ratio  of  air  to  fuel  was 
1.18  and  the  air-fuel  ratio  was  12.78.  The  overall  equivalence  ratio  was  approximately 
1.12.  The  air  and  fuel  were  injected  into  the  combustor  at  different  temperatures.  The  inlet 
temperature  of  the  fuel  was  300  K  while  that  of  the  air  was  at  589  K.  A  comprehensive 
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description  of  the  experimental  rig  and  an  analysis  of  experimental  results  can  be  found  in 
Lewis  and  Smoot,  1981. 

The  numerical  simulation  of  the  turbulent  diffusion  flame  utilized  two  different  grids,  a 
40  x  20  grid  and  a  60  x  30  grid.  A  representative  40  x  10  grid  is  given  is  Figure  32.  The 
grid  was  concentrated  in  the  region  between  the  fuel  and  air  streams  so  as  to  obtain  some 
measure  of  resolution  of  the  shear  layer  between  these  two  flowing  streams.  Such  a 
resolution  resulted  in  rather  large  cell  aspect  ratios  in  the  downstream  portion  of  the 
combustor  (Ax/Ay  =  100).  The  large  aspect  ratio  did  not  however  pose  a  problem  in 
obtaining  reasonable  convergence  for  this  flow. 

The  inlet  air  and  fuel  streams  were  assumed  to  have  uniform  velocity,  turbulent  kinetic 
energy,  and  turbulent  dissipation  rate  profiles.  A  turbulent  intensity  of  six  percent  was 
assumed  for  both  streams.  The  turbulent  length  scale  for  the  fuel  was  given  as  five  percent 
of  the  jet  diameter.  The  turbulent  length  scale  for  the  air  stream  was  assumed  to  be  ten 
percent  of  the  annular  separation  distance  (i.e.  0.1  (r0  -  rj) ).  The  numerical  solution  of  the 
flow  was  assumed  converged  when  the  mass  residual  was  reduced  4  orders  of  magnitude. 

The  comparisons  of  experimental  and  calculated  temperature  profiles  at  various  axial 
locations  (given  by  x/D  where  D  is  the  diameter  of  the  combustor)  are  given  in  Figures  33 
to  35.  The  predicted  temperatures  are  consistently  larger  than  those  given  by  experiment. 
The  60  x  30  grid  is  seen  to  yield  better  results.  For  x/D  of  0.47,  the  temperature  profiles 
are  discontinuous.  Sharp  temperature  gradients  are  possible  for  laminar  flames  but  are  not 
observed  in  the  mean  flowfield  of  a  turbulent  flow.  Turbulent  fluctuations  act  to  diffuse 
such  features  for  turbulent  flows.  This  rather  curious  behavior  is  the  result  of  the  double 
delta  function  pdf.  Jones  and  Whitelaw,  1982  found  that  a  double  delta  function  pdf  can 
yield  discontinuities  in  the  weighted  temperature  field.  Similar  observations  have  been 
made  by  other  researchers  studying  boundary  layer  diffusion  flames  (Chen,  et.  al.,  1991). 
The  calculated  temperature  profiles  indicate  that  the  temperature  increases  nearer  the  walls 
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and  with  distance  along  the  combustion  chamber.  The  experimental  observations  indicate  a 
similar  behavior  except  at  x/D  of  2.34.  The  temperature  data  at  this  location  were  not 
corrected  for  conduction  and  radiation  losses.  It  is  not  clear  how  influential  this  is  in 
making  comparisons  of  experiment  and  calculation. 

The  calculated  profiles  of  mean  mixture  fraction  given  in  Figures  36  through  39  give 
reasonable  agreement  with  experiment.  The  mixture  fraction  decreases  with  increase  in 
radius  and  with  distance  along  the  duct  This  is  consistent  with  the  increase  in  temperature 
in  these  regions.  The  form  of  the  pdf  does  not  lead  to  discontinuities  in  the  mean  mixture 
fraction.  This  was  also  observed  by  Chen,  et.  al.,  1991  for  a  turbulent  boundary  layer 
diffusion  flame. 

The  convergence  rates  corresponding  to  the  calculations  performed  in  this  simulation 
are  given  in  Figure  40.  The  60  x  30  grid  calculation  exhibits  some  initial  oscillations  but  is 
still  convergent.  The  oscillations  may  be  the  result  of  high  under  relaxation  (0.7  for  all 
variables  except  pressure  which  was  not  relaxed).  The  calculations  presented  did  not  use  a 
harmonic  mean  for  interpolating  the  exchange  coefficient.  Use  of  a  harmonic  mean  was 
detrimental  to  convergence  rate  (1200  iterations  to  converge  the  60  x  30  grid)  and  was  not 
observed  to  improve  the  calculation  in  this  case. 

5.5.  Turbulent  Premixed  Flame 

The  model  of  turbulent  premixed  flames  was  assessed  by  comparisons  with  data 
obtained  from  the  flow  in  a  sudden  expansion  dump  combustor  of  area  ratio  equal  to  2.25. 
A  propane/air  mixture  of  equivalence  ratio  0.65  was  fed  into  an  axisymmetric  sudden 
expansion.  The  inlet  Reynolds  number  was  approximately  128,000  and  the  flow  Mach 
number  was  quite  low  at  0.05.  The  inlet  flow  to  the  combustion  chamber  was 
experimentally  observed  to  be  axisymmetric. 
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The  flow  was  simulated  on  two  different  grids,  a  50  x  20  and  a  60  x  30  grid.  The  grid 
structure  was  refined  in  the  near  wall  regions  so  as  to  better  resolve  the  large  flow 
gradients.  The  inlet  velocity  and  kinetic  energy  were  specified  according  to  the 
experimental  data.  Since  the  azimuthal  fluctuations  could  not  be  obtained,  it  was  assumed 
that  the  turbulence  at  the  inlet  was  isotropic.  Therefore, 

w’w’  =  1/2  (  u'u'  +  v'v’  )  (5.1) 

The  turbulent  kinetic  energy  at  the  inlet  was  thus  determined  as  follows: 

k  =  3/4  (  u’u'  +  v'v'  )  (5. 2) 


The  turbulence  dissipation  rate,  e  at  the  inlet  was  determined  from  the  following  relation: 


Cuk3/2 

Lt 


(5.3) 


The  distribution  of  the  turbulent  length  scale  was  obtained  from  Schlichting,  1979  for  a 
fully  developed  pipe  flow: 


LT  =  R  [  0. 14  -  0.08  ( 1  -  £  )2  -  0.06  ( 1  -  £  )4  ]  (5. 4) 

The  turbulent  length  scale  was  so  specified  so  as  to  minimize  arbitrariness  in  the 
specification  of  the  boundary  conditions  for  the  turbulence  quantities.  All  calculations  of 
the  turbulent  premixed  flame  were  presumed  to  be  converged  when  the  mass  residual  was 
reduced  by  five  orders  of  magnitude. 

The  results  for  the  mean  velocity  profiles  imply  that  the  core  region  is  predicted  to  be 
essentially  inviscid  (Figures  41  to  44,  where  H  is  the  step  height).  The  uniform  velocity 
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distribution  corresponds  to  a  predicted  uniform  distribution  of  the  turbulent  kinetic  energy 
(Figures  46  to  49).  Apparently,  the  production  of  turbulence  kinetic  energy  from  the  shear 
layer  has  not  been  given  sufficient  opportunity  to  diffuse  radially  for  the  axial  locations  near 
the  inlet  (i.e.,  x/H  =  2,  x/H  =  8).  With  distance  along  the  combustor  axis,  the  mean 
velocity  is  underpredicted  (Figures  43  and  44).  The  velocity  distribution  at  x/H  equal  to  18 
seems  to  imply  that  the  combustion  is  occurring  primarily  at  the  wall.  This  is  clearly  seen 
in  the  temperature  contours  given  in  Figure  45.  The  higher  temperatures  at  the  wall  lower 
the  fluid  density  and  thus  the  velocity  increases  to  satisfy  the  continuity  constraint.  The 
experimental  results  clearly  indicate  a  greater  velocity  at  the  wall  than  at  the  centerline 
(Figure  44). 

The  profiles  of  the  turbulent  kinetic  energy  (Figures  46  through  49)  reveal  the  region  of 
high  turbulence  production  caused  by  the  shear  layer  between  the  core  flow  and 
recirculating  flow.  The  turbulent  kinetic  energy  as  given  in  Figure  46  clearly  reveals  the 
region  of  maximum  turbulence  intensity  as  obtained  from  prediction.  The  turbulence 
kinetic  energy  is  predicted  to  dissipate  with  length  along  the  combustor  whereas  the 
experiment  reveals  a  substantial  increase  in  turbulence  kinetic  energy.  This  increase  is  seen 
to  be  extremely  large  at  x/H  equal  to  18  (Figure  49).  The  region  in  which  k  increases  the 
most  is  near  the  combustor  wall.  It  is  also  within  this  region  where  most  of  the  combustion 
seems  to  be  occurring.  The  underprediction  in  the  turbulent  kinetic  energy  may  therefore 
stem  from  the  assumption  of  constant  density  flow  which  was  invoked  in  the  development 
of  the  turbulence  model. 

The  convergence  histories  given  in  Figure  50  show  that  the  calculation  of  turbulent 
reacting  flow  is  convergent.  The  calculation  on  a  50  x  20  grid  used  a  relaxation  of  0.7  on 
all  variables  except  pressure.  Pressure  was  not  relaxed  for  this  case.  The  calculation  on  a 
60  x  30  grid  also  used  a  relaxation  of  0.7  on  all  variables  except  the  pressure.  For  the  60  x 
30  grid  calculation,  the  pressure  relaxation  factor  was  0.8.  When  compared  to  the 
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convergence  histories  for  the  isothermal  turbulent  sudden  expansion  flow  (Figure  22),  it  is 
obvious  that  a  variable  density  flow  requires  much  more  computational  effort 

5.6.  Turbulent  Swirl  Flow  in  an  Axisymmetric  Sudden  Expansion 

The  results  of  a  turbulent  swirl  flow  prediction  are  very  sensitive  to  the  prescribed  inlet 
conditions  (Nallasamy  and  Chen,  1985).  Without  proper  knowledge  of  all  variables  at  the 
inlet,  predictions  may  differ  substantially  from  the  experimental  results  (Abujelala  and 
Lilley,  1983).  The  present  calculation  has  been  made  in  order  to  demonstrate  that  the 
solution  procedure  is  convergent  for  calculations  made  for  turbulent  swirling  flows.  The 
accuracy  of  the  prediction  is  a  function  of  the  turbulence  model,  accuracy  of  the 
discretization,  and  the  correctness  of  the  prescribed  inlet  conditions. 

The  experimental  results  of  Nejad,  et.  al.,  1989  were  used  for  comparative  purposes. 
The  inlet  flow  Reynolds  number  was  125,000  and  the  swirl  number  was  0.3.  In  the 
experimental  setup,  the  swirler  was  positioned  upstream  of  the  dump  plane  at  a  distance 
equivalent  to  the  inlet  pipe  radius.  The  flow  field  was  probed  for  velocities  and  turbulent 
fluctuations  using  a  two-component  four-beam  laser  doppler  velocimeter.  Further  details 
regarding  the  experimental  setup  and  the  results  obtained  from  the  experiment  are  to  be 
found  in  Nejad,  et.  al.,  1989. 

The  flowfield  simulation  of  the  turbulent  swirling  flow  was  made  using  two  different 
grids:  a  50  x  20  and  a  60  x  40  grid.  The  50  x  20  grid  is  given  in  Figure  51.  Note  that  the 
grid  structure  is  non-orthogonal.  It  was  desired  to  know  what  effect  the  non-orthogonality 
would  have  on  the  convergence  rate  for  a  swirling  flow  calculation. 

The  inlet  boundary  plane  was  selected  to  be  just  downstream  of  the  swirler.  It  was 
assumed  that  the  flow  from  the  swirler  was  essentially  a  solid  body  rotation  (Abujelala  and 
Lilley,  1983).  The  radial  component  of  velocity  was  assumed  to  be  zero.  The  inlet  axial 
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velocity  was  specified  as  uniform  as  were  the  distributions  of  turbulent  kinetic  energy  and 
dissipation  rate.  The  inlet  turbulence  intensity  was  assumed  to  be  six  percent  and  the 
turbulent  length  scale  was  specified  as  equal  to  ten  percent  of  the  inlet  radius. 

Figures  52  to  56  show  the  axial  velocity  profiles  at  various  axial  locations  (H  is  the  step 
height).  The  comparison  is  fairly  good  in  the  inlet  region  (x/H  =  4).  The  flow 
development  however  is  not  well  predicted.  The  experimental  results  reveal  that  the  flow 
at  x/H  equal  to  18  is  essentially  a  solid  body  rotation  (see  Figure  56  and  61).  The 
experimentally  obtained  axial  velocity  profile  is  seen  to  be  essentially  uniform  (Figure  56) 
while  the  swirl  velocity  varies  linearly  with  radius  (Figure  61).  The  predictions  of 
turbulent  kinetic  energy  are  qualitatively  correct  near  the  inlet  (Figure  62),  but  do  not 
develop  in  the  same  manner  as  the  experimentally  obtained  values.  Near  the  combustor 
exit,  the  experiment  reveals  an  essentially  constant  turbulent  kinetic  energy  whereas  the 
prediction  varies  substantially  with  radius  (Figure  66). 

The  convergence  histories  for  the  calculations  made  on  50  x  20  and  60  x  40  grids  are 
given  in  Figure  67.  The  relaxation  factors  used  for  the  50  x  20  calculation  were:  0.7  for  u 
and  v,  1.0  for  pressure,  0.5  for  the  swirl  velocity,  and  0.8  for  turbulence  variables  k,  e, 
and  |1t-  The  relaxation  factors  used  for  the  60  x  40  calculation  were:  0.5  for  u  and  v,  0.6 
for  pressure,  0.3  for  the  swirl  velocity,  and  0.7  for  turbulence  variables  k,  e,  and  px-  The 
solution  was  seen  to  be  convergent  for  both  calculations.  As  in  other  calculations  using 
non-orthogonal  grids,  use  of  a  harmonic  mean  for  interpolation  led  to  convergence 
difficulties.  In  some  cases,  an  otherwise  convergent  solution  would  diverge  if  the 
harmonic  mean  was  used. 
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CHAPTER  6:  SUMMARY  AND  RECOMMENDATIONS 


A  numerical  algorithm  for  simulating  turbulent  reacting  flows  has  been  developed  The 
method  is  based  on  the  solution  of  the  Favre  averaged  equations  which  govern  a  turbulent 
reacting  flow.  To  obtain  closure  of  the  governing  equations,  the  algorithm  employs  the 
standard  k-e  model  with  wall  functions  (Launder  and  Spalding,  1974).  Combustion 
models  appropriate  for  diffusion  and  premixed  flames  have  been  incorporated  into  the 
algorithm.  The  model  for  diffusion  flames  assumes  mixing  limited  combustion  and  uses  a 
double-delta  function  pdf  for  the  determination  of  the  mass  fraction,  temperature,  and 
density  distributions.  The  model  for  premixed  flames  also  assumes  a  mixing  limited 
combustion  process.  An  eddy  break  up  model  of  the  form  proposed  by  Magnussen  and 
Hjertager,  1976  is  used  to  relate  the  average  reaction  rate  to  the  mean  flow  variables. 

The  governing  equations  are  solved  numerically  by  using  a  finite  volume  method  based 
on  a  non-staggered  grid  for  velocities  and  pressure  (Rodi,  et.  al.,  1987).  The  use  of  a  non- 
staggered  grid  eliminates  the  use  of  grid  sensitive  curvature  terms.  The  governing 
equations  are  written  in  terms  of  non-orthogonal  curvilinear  coordinates  (i.e.  body  fitted 
coordinates)  and  are  discretized  using  the  hybrid  differencing  scheme.  A  decoupled  or 
segregated  solution  procedure  has  been  adopted  for  the  solution  of  the  set  of  equations 
which  govern  a  turbulent  reacting  flow. 

Tests  of  the  flowfield  simulation  technique  reveal  that  the  major  shortcomings  are  in  the 
mathematical  models  which  are  used  to  close  the  governing  equations.  Observations 
concerning  the  k-e  turbulence  model  are  consistent  with  the  findings  of  other  researchers 
(Martinuzzi  and  Pollard,  1989).  Fully  developed  pipe  flow  and  to  some  extent  sudden 
expansion  flows  are  adequately  predicted  with  a  k-e  model.  Turbulent  swirl  flows  are  not 
well  described  by  the  standard  k-e  turbulence  model. 
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Simulation  of  a  turbulent  diffusion  flame  reveals  the  inadequacy  of  a  double  delta 
function  pdf.  Discontinuities  in  the  temperature  profiles  are  clearly  unphysical  for  this 
flow.  The  calculation  of  a  turbulent  premixed  flame  in  an  axisymmetric  sudden  expansion 
resulted  in  reasonable  qualitative  comparisons  with  the  experimental  mean  velocity  profiles. 
Turbulent  kinetic  energy  was  not  successfully  predicted.  Possibly  this  is  due  to 
simplifications  made  in  the  equations  for  k  and  e. 

The  solution  procedure  itself  was  found  to  be  convergent  for  all  tests  considered. 
Proper  source  terms  linearization  was  found  to  be  an  imperative  for  obtaining  convergence. 
The  use  of  a  harmonic  mean  for  interpolating  the  exchange  coefficient  was  adopted  for 
some  tests  but  was  found  to  be  detrimental  to  the  convergence  rate. 

The  pressure  correction  equation  implemented  in  the  present  algorithm  is  appropriate 
for  coordinate  systems  which  are  non-orthogonal.  Retention  of  the  non-orthogonal  terms 
in  the  pressure  correction  equation  was  found  to  improve  the  convergence  behavior  of  the 
algorithm  for  grids  exhibiting  significant  grid  skewness  (Peric,  1990).  Although  use  of 
highly  non-orthogonal  grids  is  not  advisable  due  to  loss  of  accuracy  (Fletcher,  1988b),  it  is 
unavoidable  for  some  cases.  In  such  instances,  the  use  of  the  full  pressure  correction 
equation  will  remain  convergent. 

Based  on  the  results  obtained  from  tests  of  the  numerical  algorithm,  it  is  recommended 
that  only  low  Mach  number  flows  be  simulated  with  the  present  procedure.  Issues  and 
questions  relevant  to  finite  rate  effects  in  combustion  systems  cannot  be  addressed  with  the 
present  numerical  algorithm.  Further  improvements  of  the  present  procedure  will  involve 
increasing  the  accuracy  of  the  discretization  procedure  and  implementing  a  multi-grid  based 
algorithm.  In  the  tests  of  the  present  method,  clear  evaluations  of  the  models  used  could 
not  be  made  without  considering  grid  dependencies.  By  incorporating  a  high  order 
accurate  scheme  and  using  a  multi-grid  method  for  fine  grid  calculations,  it  is  hoped  that  the 
grid  dependency  issue  can  be  resolved.  In  this  way,  an  efficient  and  accurate  numerical 
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method  will  be  available  for  testing  the  strengths  and  weaknesses  of  various  physical 
models  and  for  simulating  the  flow  and  combustion  processes  in  modem  combustion 
chambers. 
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Table  1:  Values  of  diffusive  exchange  coefficient  and  source  term  for  variables  used  in 
simulation  of  turbulent  reacting  flow. 
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(a)  Transformed  Plane 


(b)  Physical  Plane 


Figure  1:  The  grid  in  the  transformed  plane  (a)  appears  as  non-orthogonal  in  the 
physical  plane  (b). 
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Figure  2:  The  discretization  of  the  governing  flow  equations  is  accomplished  by 
dividing  the  domain  into  finite  control  volumes.  A  sudden  expansion 
geometry  with  a  10  x  3  grid  is  given  as  an  example. 


Figure  3:  A  finite  control  volume  with  its  nearest  neighbors  is  shown.  Integration 
of  the  governing  equations  over  this  finite  control  volume  results  in  a  set 
of  non-linear  algebraic  equations. 


Calculate  cell  face  values  of  exchange  coefficient 
Calculate  coefficients  ( Ae,  Aw,  etc. ) 

Solve  for  the  u  component  of  velocity: 

1)  Determine  source  term  for  u. 

2)  Submit  to  adi  routine  for  solution. 


Solve  for  the  v  component  of  velocity: 

1)  Determine  source  term  for  v. 

2)  Submit  to  adi  routine  for  solution. 

Calculate  cell  face  fluxes. 


Solve  for  the  pressure  correction: 

1)  Calculate  coefficients. 

2)  Calculate  source  terms  for  P. 

3)  Submit  to  adi  routine  for  solution. 


Correct  u,  v,  U,  V,  and  P 


Apply  integral  mass  adjustment  if  desired 


Solve  for  scalar  variables: 

1)  Calculate  cell  face  values  of  exchange  coefficient 

2)  Calculate  coefficients  of  discretized  equation 

3)  Calculate  source  term 

4)  Submit  to  adi  routine  for  solution 


Update  thermodynamic  and  transport  properties 


Figure  4:  Order  of  events  in  obtaining  a  solution  to  the  discretized  equations 
governing  the  flow. 
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Figure  5:  Comparison  of  experimentally  and  numerically  obtained  profiles 
of  mean  velocity  at  x/D  equal  to  2  in  a  developing  turbulent  pipe  flow. 
Experimental  results  are  those  obtained  by  Richman  and  Azad,  1973. 
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Figure  6:  Comparison  of  experimentally  and  numerically  obtained  profiles 
of  mean  velocity  at  x/D  equal  to  5  in  a  developing  turbulent  pipe  flow. 
Experimental  results  are  those  by  Richman  and  Azad,  1973. 
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Figure  7:  Comparison  of  experimentally  and  numerically  obtained  profiles 
of  mean  velocity  at  x/D  equal  to  10  in  a  developing  turbulent  pipe  flow. 
Experimental  results  are  those  obtained  by  Richman  and  Azad,  1973. 
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Figure  8:  Comparison  of  experimentally  and  numerically  obtained  profiles 
of  mean  velocity  at  x/D  equal  to  15  in  a  developing  turbulent  pipe  flow. 
Experimental  results  are  those  obtained  by  Richman  and  Azad,  1973. 


89 


Figure  9:  Comparison  of  experimentally  numerically  obtained  profiles 
of  mean  velocity  at  x/D  equal  to  20  in  a  developing  turbulent  pipe  flow. 
Experimental  results  are  those  obtained  by  Richman  and  Azad,  1973. 
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Figure  10:  Comparison  of  experimentally  and  numerically  obtained  values  of 
the  mean  streamwise  velocity  component  along  the  pipe  centerline.  Experimental 
results  are  those  obtained  by  Richman  and  Azad,  1973. 


90 


Re  =  300,000 
80  x  40  grid 
r/R  =  0.6 
°  Experiment 
- Calculation 


Figure  1 1:  Comparison  of  experimentally  and  numerically  obtained  values  of 
the  mean  streamwise  velocity  component  at  a  radius  ratio  of  0.6.  Experimental 
results  are  those  obtained  by  Richman  and  Azad,  1973. 


Re  =  300,000 
80  x  40  grid 
r/R  =  0.85 
°  Experimem 
- Calculator 


Figure  12:  Comparison  of  experimentally  and  numerically  obtained  values  of 
the  mean  streamwise  velocity  component  at  a  radius  ratio  of  0.85. 
Experimental  results  are  those  obtained  by  Richman  and  Azad,  1973. 
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Figure  13:  Comparison  of  experimental  and  numerically  obtained  profiles  of 
dimensionless  turbulent  kinetic  energy  for  fully  developed  turbulent  pipe  flow. 
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Figure  14:  Comparison  of  experimentally  and  numerically  obtained 
velocity  profiles  at  x/H  of  2.25  for  turbulent  flow  in  an  axisymmetric 
sudden  expansion.  Experimental  results  from  Craig  (Craig,  et.  al.,  1984). 
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Figure  15:  Comparison  of  experimentally  and  numerically  obtained 
velocity  profiles  at  x/H  of  5.6  for  turbulent  flow  in  and  axisymmetric 
sudden  expansion.  Experimental  results  from  Craig  (Craig,  et.  al.,  1984). 
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Figure  16:  Comparison  of  experimentally  and  numerically  obtained 
velocity  profiles  at  x/H  of  7.8  for  turbulent  flow  in  an  axisymmetric 
sudden  expansion.  Experimental  results  from  Craig  (Craig,  et.  al.,  1984). 
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Figure  17:  Comparison  of  experimentally  and  numerically  obtained  velocity 
profiles  at  x/H  of  13.5  for  turbulent  flow  in  an  axisymmetric  sudden 
expansion.  Experimental  results  from  Craig  (Craig,  et.  al.,  1984). 
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Figure  18:  Comparison  of  experimentally  and  numerically  obtained  profiles  of 
turbulent  kinetic  energy  ratio  at  x/H  of  2.25  for  flow  in  an  axisymmetric 
sudden  expansion.  Experimental  results  from  Craig  (Craig,  et.  al.,  1984). 
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Figure  19:  Comparison  of  experimentally  and  numerically  obtained  profiles  of 
turbulent  kinetic  energy  at  x/H  of  5.6  for  flow  in  an  axisymmetric  sudden 
expansion.  Experimental  results  from  Craig  (Craig,  et.  al.,  1984). 
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Figure  22:  Convergence  histories  corresponding  to  calculations  of  turbulent 
flow  in  an  axisymmetric  sudden  expansion. 
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Figure  23:  Comparison  of  experimental  and  calculated  temperature  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et.  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.24. 
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Figure  24:  Comparison  of  experimental  and  calculated  temperature  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et.  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.47. 
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Figure  25:  Comparison  of  experimental  and  calculated  temperature  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et.  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.98. 
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Figure  26:  Comparison  of  experimental  and  calculated  velocity  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.24. 
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Figure  27:  Comparison  of  experimental  and  calculated  velocity  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et.  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.47. 
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Figure  28:  Comparison  of  experimental  and  calculated  velocity  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et.  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.98. 
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Figure  29:  Comparison  of  experimental  and  calculated  fuel  fraction  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et.  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.24. 
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Figure  30:  Comparison  of  experimental  and  calculated  fuel  fraction  profiles 
for  the  laminar  diffusion  flame  of  Mitchell,  et.  al.,  1980.  The  results 
correspond  to  an  axial  location  of  x/D  equal  to  0.47. 
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Figure  32:  A  representative  grid  used  for  the  calculation  of  the  turbulent  diffusion  flame  of  Lewis 
and  Smoot,  1981.  The  left  boundary  is  the  fuel/air  inlet.  The  40  x  10  grid  is  refined  in  the 
separation  region  between  fuel  and  air  thus  leading  to  high  aspect  ratios  near  the  combustor  exit. 


Figure  33:  Comparison  of  experimental  and  calculated  temperature  profiles 
for  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981.  The  profiles 
correspond  to  an  axial  location  of  x/D  equal  to  0.47. 
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Figure  34:  Comparison  of  experimental  and  calculated  temperature  profiles 
for  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981.  The  profiles 
correspond  to  an  axial  location  of  x/D  equal  to  1.96. 


104 


Figure  35:  Comparison  of  experimental  and  calculated  temperature  profiles 
for  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981.  The  profiles 
correspond  to  an  axial  location  of  x/D  equal  to  2.34. 
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Figure  36:  Comparison  of  experimental  and  calculated  mixture  fraction  profiles 
for  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981.  The  profiles 
correspond  to  an  axial  location  of  x/D  equal  to  0.47. 
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Figure  37:  Comparison  of  experimental  and  calculated  mixture  fraction  profiles 
for  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981.  The  profiles 
correspond  to  an  axial  location  of  x/D  equal  to  0.86. 
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Figure  38:  Comparison  of  experimental  and  calculated  profiles  of  mixture  fraction 
for  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981.  The  profiles 
correspond  to  an  axial  location  of  x/D  equal  to  1.96. 
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Figure  39:  Comparison  of  experimental  and  calculated  mixture  fraction  profiles 
for  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981.  The  profiles 
correspond  to  an  axial  location  of  x/D  equal  to  2.34. 
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Figure  40:  The  convergence  histories  for  the  two  calculations  made 
of  the  turbulent  diffusion  flame  of  Lewis  and  Smoot,  1981. 
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Figure  41:  Comparison  of  experimentally  and  numerically  obtained 
axial  velocity  profiles  at  x/H  equal  to  2  for  a  turbulent  premixed  flame. 
Experimental  results  from  Ahmed  and  Nejad,  1990. 
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Figure  42:  Comparison  of  experimentally  and  numerically  obtained 
axial  velocity  profiles  at  x/H  equal  to  8  for  a  turbulent  premixed  flame. 
Experimental  results  from  Ahmed  and  Nejad,  1990. 
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Figure  43:  Comparison  of  experimentally  and  numerically  obtained 
axial  velocity  profiles  at  x/H  equal  to  12  for  a  turbulent  premixed  flame. 
Experimental  results  from  Ahmed  and  Nejad,  1990. 
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Figure  44:  Comparison  of  experimentally  and  numerically  obtained 
axial  velocity  profiles  at  x/H  equal  to  18  for  a  turbulent  premixed  flame. 
Experimental  results  from  Ahmed  and  Nejad,  1990. 
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Figure  46:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  at  x/H  equal  to  2.  Experimental 
results  from  Ahmed  and  Nejad,  1990. 
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Figure  47:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  at  x/H  equal  to  8.  Experimental 
results  from  Ahmed  and  Nejad,  1990. 


Ill 


r  /  R 


°  Experiment 

-  50  x  20  grid 

-  60  x  30  grid 


Figure  48:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  at  x/H  equal  to  12.  Experimental 
results  from  Ahmed  and  Nejad,  1990. 
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Figure  49:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  at  x/H  equal  to  18.  Experimental 
results  from  Ahmed  and  Nejad,  1990. 
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Figure  50:  Convergence  histories  for  the  calculations  made  of  the 
turbulent  premixed  flame  of  Ahmed  and  Nejad,  1990. 
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Figure  51:  Grid  distribution  used  in  calculating  turbulent  flow  in  an  axisymmetric  sudden 
expansion  with  and  without  swirl. 
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Figure  52:  Comparison  of  the  calculated  and  experimental  axial  velocity 
profiles  for  a  turbulent  swirling  flow.  The  results  are  presented 
for  x/H  equal  to  4.  Experimental  results  from  Nejad,  et.  al.,  1989. 
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Figure  53:  Comparison  of  experimentally  and  numerically  obtained  axial 
velocity  profiles  for  a  turbulent  swirling  flow.  The  profiles  correspond 
to  x/H  equal  to  6.  Experimental  results  from  Nejad,  et.  al.,  1989. 
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Figure  54:  Comparison  of  experimentally  and  numerically  obtained 
axial  velocity  profiles  for  a  turbulent  swirling  flow.  The  profiles 
correspond  to  an  axial  location  of  x/H  equal  to  8.  Experimental 
results  from  Nejad,  et.  al.,  1989. 
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Figure  55:  Comparison  of  experimentally  and  numerically  obtained 
axial  velocity  profiles  for  a  turbulent  swirling  flow.  The  profiles 
correspond  to  an  axial  location  of  x/H  equal  to  10.  Experimental 
results  from  Nejad,  et.  al.,  1989. 
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Figure  56:  Comparison  of  experimentally  and  numerically  obtained 
axial  velocity  profiles  for  a  turbulent  swirling  flow.  The  profiles 
correspond  to  an  axial  location  of  x/H  equal  to  18.  Experimental 
results  from  Nejad,  et.  al.,  1989. 
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Figure  57:  Comparison  of  experimentally  and  numerically  obtained 
swirl  velocity  profiles  for  a  turbulent  swirling  flow.  The  profiles 
correspond  to  an  axial  location  of  x/H  equal  to  4.  Experimental 
results  from  Nejad,  et.  al.,  1989. 
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Figure  58:  Comparison  of  experimentally  and  numerically  obtained 
profiles  of  swirl  velocity  for  a  turbulent  swirling  flow.  The  profiles 
correspond  to  an  axial  location  of  x/H  equal  to  6.  Experimental 
results  from  Nejad,  et.  al.,  1989. 
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Figure  59:  Comparison  of  experimentally  and  numerically  obtained 
swirl  velocity  profiles  in  a  turbulent  swirling  flow.  The  profiles 
correspond  to  a  location  of  x/H  equal  to  8.  Experimental  results  from 
Nejad,  et.  al.,  1989. 
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Figure  60:  Comparison  of  experimentally  and  numerically  obtained 
swirl  velocity  profiles  for  a  turbulent  swirling  flow.  The  profiles 
correspond  to  an  axial  location  of  x/H  equal  to  10.  Experimental  results 
from  Nejad,  et.  al.,  1989. 
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Figure  61:  Comparison  of  experimentally  and  numerically  obtained 
swirl  velocity  profiles  for  a  turbulent  swirling  flow.  The  profiles 
correspond  to  an  axial  location  of  x/H  equal  to  18.  Experimental 
results  from  Nejad,  et.  al.,  1989. 
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Figure  62:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  for  a  turbulent  swirling  flow.  The 
profiles  correspond  to  x/H  equal  to  4.  Experimental  results  from 
Nejad,  et.  al.,  1989. 
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Figure  63:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  for  a  turbulent  swirling  flow.  The 
profiles  correspond  to  x/H  equal  to  6.  Experimental  results  from 
Nejad,  et.  al.,  1989. 
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Figure  64:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  for  a  turbulent  swirling  flow.  The 
profiles  correspond  to  an  axial  location  of  x/H  equal  to  8.  Experimental 
results  from  Nejad,  et.  al.,  1989. 


Re  =  122,728 
S  =  0.3 

°  Experiment 
50  x  20  grid 
60  x  40  grid 


Figure  65:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  for  a  turbulent  swirling  flow. 

The  profiles  correspond  to  an  axial  location  of  x/H  equal  to  10. 
Experimental  results  from  Nejad,  et.  al.,  1989. 
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Figure  66:  Comparison  of  experimentally  and  numerically  obtained 
turbulent  kinetic  energy  profiles  for  a  turbulent  swirling  flow. 

The  profiles  correspond  to  an  axial  location  of  x/H  equal  to  18. 
Experimental  results  from  Nejad,  et.  al.,  1989. 


Figure  67:  Convergence  histories  for  the  two  calculations  made  of 
a  turbulent  swirling  flow  in  an  axisymmetric  sudden  expansion. 
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APPENDIX  A 


The  equations  which  are  taken  to  govern  laminar  flow  in  axisymmetric  geometries  (r,z)  are 
given: 
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,r2v  2  /I  d(rv)  ,  d  uV 
,89=i‘[t-3(i  dr  +  dz] 

(A.  6) 

*>|  n' 
as  1^° 

+ 

«SI* 

1 _ 1 

=1 

II 

*P 

(A.  7) 

S  (Puh°> +  F5  (PIYh«)  ■  k  fe  Tte] + [r  ft  T5?] '  h° 

1^1 

3  ,  v  >  13,  a  I"  dYil  1  3  r  H  3  Ysl  . 
di  (puYi) + 7a?  (prvYi) = LsH  -arJ+^Lr  sS  arJ  +  toi 


(A.  8) 

(A.  9) 
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APPENDIX  B 


The  equations  which  are  taken  to  govern  turbulent  flow  in  axisymmetric  geometries 
(r,z)  when  an  eddy  viscosity  model  is  used  to  close  the  system  of  equations  are  given: 


|(pu)  +  i|(Prv)  =  0 


(B.l) 


^  (PSO)  +  (prvu)  =  -  H  ^  [reff  ^  [rT*f{  I?]  (B-2) 

|(pnv)+||(pr^)  =  -|  +  |[reff|]  +  l|[rreff|]  (B.3) 


^  (p®y + tJ?  (pr^} 


ii 

r3r 


rTeff  9h7 
Prx  dr 


N 

■1-S 

i=l 


© 


»? 


(B.4) 


+  H(P^i)=| 


O  _  I  _  'V_  O  _  *V  T  U 


Scx  dz  J  +  r3r[Sc7  dr  J 


©i 


(B.5) 


|  (PBk)  +  TS  (P**>  =  |  [reff  |]  +  H  [rFeff  |]  +  Pk  -  P*  (B.6) 

|  <P5e)  + (ptfe)  =  |  [r=«  H  +  7S  [rr«  I]  + 1  <  Clpk  •  W  <B7> 

I  <pa?)  +  H  (P^>  =  |[retf  g]  +  ||[rretf  g]  +  S4(ZJ)  (B.8) 


APPENDIX  C 


The  general  conservation  equation  in  (r,z)  coordinates  is  transformed  to  general 
curvilinear  coordinates.  The  transformation  results  in  the  following  governing  equation 
and  transformation  metrics: 


|(pu«+^(pv«=J[rqn|]+|.[rq22^ 

+Mrq,2S+yr,2iS 

+  r  1 J 1  S4»  (4,ti) 

(C.  1) 

U  =  r(ur11- v^) 

(C.  2) 

V  =  r  ( vz^  -  ur^ ) 

(C.  3) 

qn  =  i  j  |(rn2  +  zn2) 

(C.  4) 

T22  =  |  j  |  (  r£2  +  z£2  ) 

(C.5) 

<121  “  ^12  ~  *  |  j  |  (r£rti  +  ) 

(C.  6) 

1  J  1  —  (  Z^  -  Z^rj;  ) 

(C.7) 

125 


